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ABSTRACT 

The goal of the investigation was the development of a class of 
control systems which efficiently results m high-accuracy (:£ 10 radiana) 
earth -pointing attitude motions of satellites of different configurations 
m elliptic orbits. The earth-pointing orientation, 1 e , the orientation 
such that one body-fixed axis is parallel to the local vertical and another 
is normal to the orbit plane;, is required for the lifetimes of the satel- 
lites Gas jets provide the control torque 

Linear differential equations with time -varying coefficients., which 
include terms for the gravity torque due to an oblate earth and terms 
for the aerodynamic torque, are used to describe the attitude motion when 
a satellite is practically earth -pointing Nonlinear equations with time- 

varying coefficients are used to describe the attitude motion when acquisi- 
tion of the earth-pointing motion from large deviation angles (w80°) is 
considered. 

Pontryagm's Maximum Principle, the necessary conditions for exact 
solutions of optimal bounded-phase -coordinate problems and guidelines 
obtained from the minimum-fuel station-keeping controls devised for 
single-axis systems are used m the development of the station-keeping 
part of the control system The acquisition part of the control system, 
developed here, results m acquisition of the earth-pointing motion from 
large angles m the time of one-quarter orbit with comparatively little 
fuel expenditures 

The motions of the satellites with the developed station-keeping 
control systems are simulated on an analog computer, and, the performances 
of the systems are evaluated The nonlinear differential equations which 
include the developed acquisition control systems with time delays are 
integrated by using a digital computer The fuel expenditures and the 
times of acquisition obtained from these digital computer runs are com- 
pared with those for the satellite motions described approximately by 
linear differential equations with the optimal controls obtained from the 
maximum principle 


m 



The simulation studies together with the performance evaluations 
showed that the control systems are quite efficient as well as reliable 
(The acquisition system is a back-up system to the station-keeping 
system ) The overall control system, which is very simple to realise, 
is given m diagram form 
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I INTRODUCTION 


A. STATEMENT OF TEE PROBLEM 

At this time there are many planned artifical earth satellites which 
must perform tasks which require a given satellite -fixed line to he 
(nearly) parallel to the local vertical while the angles between every 
satellite -fixed line and tne orbit plane are (nearly) constant. This 
required orientation is called (near) earth-pointing Earth-pointing is 
not a natural satellite orientation so that control effort must be spent 
to keep the orientation (called station-keeping) once it has been 
acquired. The time interval of station -keeping for practical satellites 
ranges from minutes to years. In any case the controller, which causes 
the satellite to acquire the earth pointing attitude and keep it, must 
be efficient as well as reliable 

The problem considered m this investigation is how to efficiently 
and reliably control a satellite's attitude motion so that it is very 
accurately earth-pointing for a given interval of time 

B SATELLITE CONFIGURATIONS 

Since attitude motion characteristics vary widely with satellite 
configuration and since the efficiency of a given control varies with 
the characteristics of attitude motion, it, at first, seems necessary 
to consider an almost infinite variety of satellite configurations 
However, upon further investigation ic is found that a few satellite 
configurations exhibit a wide variety of attitude motions whose charac- 
teristics are basic characteristics of a whole range of satellite configurations. 

Four configurations of satellites and tneir orbits are considered 
here These four configurations were chosen since they exhibit a wide 
variety of natural attitude motions and since they are similar to some 
future as well as some recent earth satellites 

Satellites (l) and (4) are similar to some of today's "stable" 
scientific satellites e g , satellites of the Orbiting Geophysical 
Observatory and the Tiros series Satellite (4) is such that aerodynamic 
forces significantly affect its attitude motion Satellite (2) is 
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similar to "unstable" manned and unmanned spacecraft which will be 
inserted into (nearly) circular earth parking orbits for transfer to 
orbits about other celestial bodies. Examples of such spacecraft are 
Apollo and Unmanned Mars Excursion Vehicle Satellite (3) is similar to 
future "unstable" military applications satellites 

The orbits of the satellites are elliptic, although for satellites 
(2) and (3) the eccentricity is assumed to be only 0 01. The orbits 
are not considered circular since sucn an assumption is generally a 
gross oversimplification and since the number of orbits a satellite will 
complete before the aerodynamic forces cause it to move on a trajectory 
back to earth increases with eccentricity (Breakwell and Koehler [4] 
have derived an expression for the number of oroits before decay This 
expression, which gives s lifetime of about 10,000 orbits (about 1 7 years) 
for an initial perigee neight; of 200 miles and an initial eccentricity 
of 0.03, has been experimentally verified ) 

C CONTROLLERS AND PERFORMANCE CRITERIA 

In recent years many devices which showed promise of controlling the 
attitude motion of satellites have been constructed, analyzed and to 
some extent applied with success These controllers are categorized 
according to tneir subsystems which supply the control torques Control 
torques are obtained by devices such as gas jets, ion propulsion units, 
three -axes-gyroscopes, v-gyroscopes, reaction wheels, reaction spheres, 
extendable Dooms and magnetic devices 

Several investigators nave considered such devices m solving various 
acquisition and station -keeping problems For example, Eorwitz [19] has 
considered pulse width modulation for fairly high accuracy control of 
highly idealized attitude motion. Euston [20] considered the feasibility 
of employing twin-gyroscopes for control, and, Schwartz [32] considered 
minimum energy acquisition using reaction wheel control Asymptotic 
stabilization of the spin axis of a satellite m a circular orbit by 
magnetic control was considered by Wheeler [34] Eaefner [15] and 
Nichoi [28] each considered some of the general aspects of attitude 
control Although these solutions contribute to the solution of the 
problem at hand, none result m sufficient guidelines for the construction 
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of a device which will efficiently and reliably control the attitude 
motion for high accuracy earth-pointing 

Because of their simplicity and reliability, gas jets are considered 
m this investigation as the suppliers of control torque 

The error m orientation for the high accuracy earth-pointing 
orientation is considered'* to have priority over all other performances 
measures. However, for efficiency, the weight of the fuel spent for 
control must he as near minimal as is practicable (See Section B, 

Chapter III ) Another criterion used m the present investigation is the 
time limit for acquisition of the high accuracy earth-pointing orientation 
Since some satellites require rapid acquisition (eg, the acquisition of 
the Gemini -Agena to earth -pointing before firing the agena rocket for 
changing orbit), the time of acquisition is limited to the time of one- 
fourth of an orbit 

D PREVIOUS CONTRIBUTIONS 

The present investigation is a continuation of the work at Stanford 
on the attitude control of satellites Busch** [ 6 ] has developed a sub- 
optimal minimum-fuel acquisition control law for a "stable" satellite m 
an elliptic orbit The use of this control law does not consistently 
result m acquisition to high accuracy earth-pointing unless the eccen- 
tricity of the orbit is very small However, Busch has included a 
reaction wheel for improving the accuracy when the eccentricity is not 
small Although the reaction wheel must be slowed periodically by apply- 
ing gas jet torque, Busch's controller performs quite well m station- 
keeping. The controller gives fairly high accuracy m earth -pointing, 
but the controller is not efficient while station-keeping 


*After completing this work it came to our attention that A E Pearson 
m his presentation on "Performance Maintainability m Precision Attitude 
Control Systems" m J A C C preprints (Philadelphia, 1967) discussed 
generalized performance measures for practical attitude control systems 

**Also, see SUDAAR No 26 l or R E Busch and I. Flugge-Lotz, "Attitude 
Control of a Satellite m an Elliptic Orbit", Journal of Spacecraft and 
Rockets, Vol 4 , No 4 , 1967. 
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Others who have investigated certain aspects of minimum-fuel 
attitude control are Flugge-lotz and Marbach [15], Foy [l4J, Meditch [ 27 ], 
Athans [1], Craig and Flugge-Lotz [10] and Hales [1 6 ] 

E. CONTRIBUTIONS 

In Chapter II the equations of motions of rigid satellites m 
elliptic orbits about an oblate body with atmosphere are derived m terms 
of the "three-axes Euler angles" which for small angles are nearly the 
yaw, roll and pitch angles Ihe effects of disturbances such as those 
due to solar radiation pressure, meteoroids, etc. are discussed A means 
of determining an upper bounds on the error due to simplifying the full 
nonlinear equations of controlled motion is given 

In Chapter III the theory of the control of satellite attitude motion 
is discussed and its application justified 

In Chapter IV a high-accuracy station -beeping feedback control law 
is developed The application of this control law to "unstable" as well 
as "stable" satellites results m very little fuel expenditure The 
reliability and efficiency of the station-keeping controller under adverse 
circumstances is investigated m Chapter VI where the overall system is 
discussed 

Suboptimal acquisition control laws for the four satellites are 
developed m Chapter V by extending Busch's results with the aid of phase - 
plane methods and Pontryagm's Maximum Principle. 

It is assumed that sensors of attitude angles and their rates for 
such high accuracy will be available m the very near future 
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II. THE EQUATIONS OF CONTROLLED SATELLITE ATTITUDE MOTION 


Satellite attitude motion is suitably described by a system of 
differential equations m -which the dependent variables are the gener- 
alized coordinates of the satellite. Three of the generalized coordinates 
completely determine the attitude of a rigid satellite m some reference 
frame. The time -history of these three coordinates and their rates for 
some interval of time is called the attitude motion m that time interval 
or more simply the motion 

A THE INERTIA TORQUE 

For a rigid* satellite, B, the inertia torque for the mass center 
of B, B*, is well known. (See, for example, Kane [22].) If n , 
i = 1,2,3, are elements of a right-handed set of mutually perpendicular 
unit vectors which are parallel to principal axes of inertia of B for B*, 
say L^, l = 1,2,3, (see Figure 2 l) and if I , l = 1,2,3, are the 
associated principal moments of inertia, then the inertia torque of B 
for B* is given by 

Si = [ < V» 3 ( X 2 - V - “ibki 

+ tm 3 * 1 (l 3 - Ij ) - i> 2 I 2 ]n 2 (2 1) 

+ ^1^2 ^1 ~ X 2^ ~ cd 3 i 3 ]H 3 , 

where cjo^, 1 = 1,2,3, are the measure numbers of the angular velocity 
of B m an inertial reference frame for the basis, n , 1 = 1,2,3, and 
where the symbol ( ) = d( )/dt 

One orbit of the satellites considered very nearly lies m a plane 
in an inertial space. ( See King-Hele [23].). This plane, which contains 
the mass center, E*, of the earth, E, (together with the plane’s 
normal) determines a suitable inertial reference frame 


Causes of non-rigidity of the satellite and their effects on the con- 
trolled motion of the satellite and on the fuel cost are discussed m 
Section C of Chapter VI. 
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ORBIT 

SEGMENT 



Figure 2 1 Principal Axes of Inertia.* L > i = 1*2* 3* the Unit 
Vectors n^* 1 = 1*2* 3* the Satellite B and the 
Earth E 



In terms of the "three-axes Euler angles", Q^, i = 1*2,3,* (see 
Busch [6]) the measure numbers of T-j., (2 l), fo r the basis 

l = 1,2,3, axe 


T I1 = { “ e i C 2 C 3 - Vs + e i 0 2 S 2 C 3 ' ®2 0 3 C 3 + ®1 0 3 C 2 S 5 
- e + S;L s 2 c 3 ) - egC^Cg + + ^s^)] 

“ ^ s i 8 5 ‘ C 1 S 2°3^ “ ^l^^ X I 1 
T I2 = {_0 2 C 3 + 9 1 C 2 S 3 " G 1 0 2 S 2 S 3 + 0 2 0 3 S 3 + 6 1 9 3 C 2 C 3 

- 9[0 1 (c 1 c 3 - S-^Sg) + + ®g( c 1 S 2 C g - s ! s 3 ) ] 


(2.2a) 


- 0 ( S 1 C 3 + C 1 S 2 S 3 ) " W X h 
T 13 = K - 0 1 S 3 - \ e 2 C 2 + 0(0 1 S 1 C 2 + 0 2 C 1 S 2 ) 

“ 0C 1 C 2 “ k 3 £ °l“2 } X X 3 

where c = cos 9^, = sm 0^, i = 1,2,3, and where 

<\ - (8=! * 0 2> s 3 ' (9 °1 B 2 - 0 1 C 2 )O 3 
a> 2 = (9=3. + 9 2 )0 3 + ( ' iC 3 3 2 ‘ e i C 2 ,S 3 
m 3 * e °l°2 + 0 1 S 2 + 6 3 


(2 2b) 


(2.2c) 


(2 3) 


The angle 0 is defined m Appendix A (see Figure A l) 
parameters k^,k 2 ,k 3 are defined by 



and the inertia 


(2.4) 


*When 0 , l = 1,2,3, aie small, they can be considered as the yaw, roll 
and pitch angles, respectively 
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B THE (ACTIVE) TORQUE 


The torque for B* is due to all contact and body forces actxng on 
B The body forces are due primarily to the gravitational attraction of 
the earth and other celestial bodies and the interaction of the satellite 
’With the earth’s magnetic field The contact forces are due primarily to 
the interaction of the satellite with the earth’s atmosphere and emissions 
from the sum, to collisions of the satellite with meteoroids and to the 
controls, for example gas jets 

All of these forces were considered m the derivation of the torque 
for B* (See Section F of this chapter,) However, the torque expression 
used m the derivation of the equations of motion is a result of consid- 
ering only aerodynamic, control and earth gravitational forces This 
torque expression can be written as 

3 3 

T = ) Tn =)(T + T +C)n (2.5) 

Zj 1-:l Zj gl 31 1 _1 

i=l i=l 

where for the basis n , l = 1,2,3, the measure numbers C , l = 1,2,3, 

are for the control torque vector, the measure numbers T , l =1,2,3, 

gi 

are for the gravitational torque vector and the measure numbers T , 
i = 1,2,3, are for the aerodynamic torque vector. 


C. THE FULL NONLINEAR DIFFERENTIAL EQUATIONS OF CONTROLLED MOTION 

If for each n , l = 1,2,3, the corresponding measure numbers of 
Tj. and T are divided by I , l = 1,2,3, respectively, and summed to 
zero, then from (2.2a)-(2 2c) and (2 5) the result is represented by 
three equations m 0 , i = 1,2,3 If these three equations are solved 
for 0 , i = 1,2,3, they give 


®1 * ( V2 S 2 ‘ ®2®3 - Sl TT'2 + 09 2 C 1°: 


- 00 3 S 1 + 80 1 S 2 + S 3 k 2“A - 


+ C 3 T lTl " ^Vl 2 )/= : 


(2 6a) 
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®2 - S i e 3°2 - 0B 1°1 - 0S 3 O 1 B 2 ' 0S 1 


- ' s 3 k r2 m 3 

+ 0 3 I 2 /I 2 - SjV 1 ! 


(2.6l3) 


°3 = - S 1 S 2 - 0 1 0 2 C 2 + 00 lV2 + 00 2 C 1 S 2 


- 0= 1 e 2 - k^ + t 5 /i 3 


(2 6c) 


Equations ( 2 . 6a) -(2 6c) axe the desired differential equations of 
motion m the dependent variables 0 , 1 =* 1,2,3, and the independent 
variable t if 0 and 9 are known functions of t and T , 1 = 1,2,3, 
are known function of 0,0, 1 =1,2, 3, and t. Except for C , 

3. jL 1 

1 = 1,2,3, which are to be determined m the following chapters as 
functions of 0^, 0^, 1 = 1,2,3, and t, the measure numbers T^, 

1 = 1,2,3, are known functions of 0 , 0 , 1 = 1,2,3, and t Equations 

(B.ll) with (B.3) and (B.lO) give the desired relationships for T , 

al 

1 = 1,2,3, and Equations (A. 8) with (A 9 ) and (A 10 ) give the desired 
relationships for T , 1 = 1,2,3, if r, r , q/r , c 0 and Sg are 
known functions of time Equations (C 23)-(c 29) give r, q/r 3 , 6, r 2 , 

0 , s 0 and Cg* as functions of time. 

The measure numbers of the torque as given by (2 5 ) can be written 
with the aid of Equations (A.8)-(A,10), (B.3), (B 10) and (B 11) as 

T 1 - ®1 + - J 2 )l - 5C 2 S 2 B 3 


+ J( r E /r;'[5(7st k - 1 ) e 2 s 2 c 3 " 2s R c f ,c„o?o„o. 


8 5 0 1 2 3 


* 

If, 


“ 2s 5Wl S 2 C 3 + 2s 5 e 0 C l S l C 2 C 3 

- 2C & C 1 S ;L C 2 C 3 *=■ SCgSgCgC^gSj 

O 

+ lOSgCgCgC^CgSg + (TERMS IN s 
of course, 0 Q is replaced with 9 


> • 


* 

P 



(continued) 
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- (C D /2)SpV{rU 2 s 2 + £ 3 c 2 s 3 ) 

- r(0 - 0 e c b)^ 2 s 1 C 2 + & 3^ C 1 C 3 " S 1 S 2 S 3^ 

+ r ®E S 5 C 0 ^Q C 1 C 3 ~ £ 5^ a l°5 + C 1 S 2 S 3^^ ^ 2 7a ^ 

T 2 = C 2 + v/r 5 ^ - I 3 ){3c 2 s 2 c 3 + J(r E /r) 2 [5(l - 7s 2 s 2 )c 2 s 2 c 3 

" 2s 8 C & C 0 C 1 C 2 S 3 + 2C 6 C 1 C 2 S 2 C 3 + ^ S 5 C & S 0 C 1 C 2 C 3 “ ® S 8 S 0 C 0 S 1 C 2 C 3 
+ 18s| S 2 c 2 s 2 c 3 + (TEEMS IN s 2 , s 3 )]} - (C D /2)Sp Vtr^c^-^) 

f r(0 - 0 E c 5 )[^( Sl s 2 c 3 + c^g) + Vl C 2 ] 

+ r0 E C 0 S S [ ^3 (s i S 3 " C 1 S 2 C 3 } - Vl C 2 ]} (2 7b) 

T 3 = C 3 + ^2 ” 3C 2 C 3 S 3 + J ( r E //r ^ t5(7s g s 0 - l) c 2 C 3 S 3 

, q 2 2,o 2 q 2 2 2 

+ 8s & s 0 c 0 c 1 c 2 c 3 + 8s g c 5 s 0 s ;L c 2 c 3 - l8s 5 s e c 2 c 3 s 3 

- 2s 2 c 2 c 2 c 3 s 3 + 2s g c s c 0 c 2 s 2 c 2 + (TERMS IN s 2 , , ) ] } 

" (C D /2)SpV{-r(i lC2 s 3 + ^c 2 c 3 ) 

+ r < 0 “ 0 E C 5 )U 1 (C 1 C 3 - S 1 S 2 S 3 } " V S 1 S 2 C 3 + C 1 S 3 )] 

+ r0 E C 0 S & U l (s i C 3 + C 1 S 2 S 3 } - *2 (b 1 8 3 “ C 1 S 2 C 3 )]} (2 7c) 

where 

P = P p exp[K(r - r p )] (2 8) 

V =(h/a)(l + 2ec 0) 1 / 2 (2 9) 

Equations (2.6a)-(2 6c) with Equations (2.3), (2 7a)-(2 9) 

and with Equations (C.23)-(C.29) are the desired nonlinear differential 
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equations of controlled motion This model of the controlled motion can 
he simplified without introducing significant error The simplifications 
strongly depend on the parameters of the satellite and on the mode 
(steady-state or acquisition) of motion 

D PARAMETERS OF THE SATELLITES AND THE SIMPLI CATION OF THE MOTION 
EQUATIONS 

In the Introduction, Chapter I, the general configurations and the 
altitudes of the four satellites to be considered were discussed Here 
this discussion is translated into numbers for I^I^I^k^k^k^e^S, 

9 , p , the height of perigee (ph) and the height of apogee (ah) These 

P p . 

numbers for the four satellites are given m Table (2 l) together with 
the air density at apogee, p o 


TABLE 2 1 PARAMETERS OF THE SATELLITES AND OF THEIR ORBITS 


Satellite No 

( 1 ) 

( 2 ) 

(3) 

00 

^1 

2 30xl0 5 

l 15 X 10 5 

7 OOxlO 3 

1 43X10^ 


2 33xl0 5 

7 OOxlO 3 

1 21X10 5 

1 44X10^ 

X 3 

2 36x10 5 

1 2 lxl 0 5 

l 15X10 5 

l 45x10^ 

(slug-ft 2 ) 





k l 

0 01 

o 99 

-0 86 

0 01 

k 2 

-0 03 

-0 86 

-0 89 

-0 01 

k 3 

0 01 

-o 89 

0.99 

0 01 

e 

0 05 

0 01 

0 01 

0.03 

a(miles) 

4480 

4500 

4500 

4260 

ph (miles) 

300 

500 

500 

200 

ah (miles) 

7^5 

585 

585 

460 

0^ (radians) 

variable from 0 to 2 ji 



6 (radians) 

variable from 0 to it if 0 

VI 

a 

VI 


0 

6 2xl0" 15 

1 6xio _l6 

1 6xlO- 16 

■m 

P a 

3 4xio“ l8 

4 7 xio" 17 

4 7xlO” 17 


(slugs/ft 3 ) 




1 
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The orbits of the satellites are fixed by the values given for 

and 5 The shapes and compositions of the satellites are not 
completely fixed by the values given for 1^, 1^ a 11 * 3 - 1^ and are even 

less fixed if values are only given for k^, i = 1,2., 3, (since common 
multiples of 1^, l = 1,2,3, result m the same values of k , l = 1,2,3) 
Topical configurations of the four satellites are as follows 

(1) Nearly spherical with a weight of 60,000 lbs , a maximum 
dimension of 14 ft and a specific weight of 40 Ibs/ft 

(2) Nearly circular cylindrical with a weight of 50,000 lbs , 
a height of 30 ft., a mean diameter of 6 ft. and a specific 
weight of 60 Ibs/ft The axis of minimum moment of inertia 
is nearly tangent to the orbital path when m the earth- 
pomting mode' 

(3) Similar to ( 2 ) except that its axis of minimum moment of 
inertia is nearly coincident with the local vertical when m 
the earth-pointing mode 

(if-) Nearly spherical with a weight of 10,000 lbs , a maximum 
dimension of 9 ft* and a specific weight of 30 lbs. /ft 

The simplifications of the equations of motion are naturally divided 
by the regimes of [0^1, 1 = 1,2,3, and the magnitudes of their com- 
patible rates In the acquisition mode |0 |, 1 = 1,2,3, vary from 

-3 1 

about one radian down to 10 radians or less with most of the time of 

acquisition spent with |0.J, 1 = 1,2,3, assuming the larger values 

In the steady- state mode the motion will be controlled m a manner such 

that (0.^1, 1 = 1,2,3, will be less than 10 3 radians In the following 

two parts of this section. Parts 1 and 2, the equations of controlled 

motion are simplified for the two modes of motion. The reasons for the 

simplifications are the solution of the full nonlinear equations of 

controlled motion with the aid of a digital computer is extremely costly, 

the analog simulation of these equations (even after they have been 

♦ . 

linearized m 0 , 0^, 1 = 1,2,3) requires a greater number of operational 
amplifiers than is available on the analog computer used (two pace TR-^8 
computers slaved together), and, finally, the error due to the simplifica- 
tions is insignificant m the final results (see Chapter VI ) 
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1 The Simplified Acquisition Equations of Motion 


For some satellites certain terms m the nonlinear equations of 
controlled motion are insignificantly small when the controlled satellite 
is m the acquisition mode. These terms are less than one-one hundredth 
of the other terms m magnitude for all but about the last one-tenth of 
the time of acquisition, which is assumed to be the time of one -quarter 
of an orbit or less During the last one -tenth of the acquisition time 

_s r 

all of the terms are of the order of 10 or less 

For satellites (l),(2) and (3) the insignificant terms are 
those m equations (2 7a) -(2 7c) with and J as coefficients and 

some of those m equations (2 6a)-(2 6c) and 7c) which are 

products of two or three of s ,9 , l = 1,2, 3, and 9 The cost of the 
digital computer solution of the equations with the latter of the above 
insignificant terms included is not significantly greater than the cost 
of the solution without these terms If only those terms with and 

J as coefficients are omitted, then the simplified acquisition equations 
of motion for satellites (l),(2) and (3) are 


9 1 ‘ (9 1 9 2 S 2 - S 2 8 3 - 00 l s l s 2 + 09 2 c 1 c 2 ‘ 99 3 S 1 + 9c l B 2 

+ Wft ■ + - s iW' ! )v 2 v 3 


S 3 C 2^2 ' '"^‘2 ( M,/r ^ e 2 S 2 C 3 S 3^ c 2 


(2 10a) 


9 2 = 0 1 0 3°2 - "l°l - 90 3°1 S 2 - 0S 1 - Vm - k l s 3 a) 2“3 
+ o 3 C 2 /l 2 + 3k 2 (|i/r 3 )c 2 s 2 o| 


' “s 0 ]/*! + 5k 1 (t*/>’ 3 )o 2 s 2 s‘ 


9 3 " - 0 1 S 2 - 9 1 0 2 c 2 + "l s l°2 + 09 2 C 1 3 2 " 9C 1°2 ‘ k 3 ffi l“2 


+ Cjl, - Sk„((i/r 3 )e 9 o 


o 3 


2 3 3 


(2 10b) 


(2 10c) 
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•where o^, x = 1,2,3, are given by (2.3) and |i/r 3 , 9 and 0 are given 
by (C 2k ), (C 25) and (C 27) with 0 q , the angle from perigee at t = 0, 
added to nt 

For satellite (4) the terms which are insignificant and will be 
omitted m the acquisition equations of motion are those which contain 
one of JjC.-.f and x8Ji s s. ,i*j, k = 1,2,3,, as factors In this 
case, if Yrd is replaced by (h 2 /a )(l + 2ec0) (h /a )[l + 2e cos 
(nt + 0 )] (see equations (B.10) and (C.29)), the simplified equations 
of motion are the same as (2 10a) -(2 10c) except that the terms 

rhs^ = (C D /2)Sp(h 2 /a 2 )[l + 2e cos(nt + 0 q )3 

* <'A c 2 )C Vl°2 + Vl°3> (2 lla) 

rhs^ = 0 , (2 lib) 


rhs 3 = (C /21 3 ) Sp(h 2 /s3[l + 2e cos(nt + 0 )] 
X ^2 C 1 S 3 " ^1 C 1 C 3 ) 


(2.11c) 


are added to the right-hand sides of equations (2 10a), (2.10b) and (2 10c), 
respectively 

If (C 23) is substituted into (B 2) and 0 is assumed to be 0 
at t = 0, then 


p = p p exp{Ke[cos(nt + 0 q )]} 


(2 12 ) 


A more convenient form of the simplified acquisition equations 
of motion is obtained by letting r = nt, ( )'= d()/dr 

(2 13) 


*2X1-1 " A 


x 0 = 0 1 
2X1 l 


for 1 = 1, 2 , 3 


With these substitutions and the substitution of (C 2k ), (C 25) and (C 27) 
into (2 10a)-(2 10c) the equations of motion for satellites (l),( 2) and 
(3) becomes 



(continued) 



*2 - (* 2 V2 - X 4 X 6 + V X 4 C 1 C 2 - X 2 S 1 S 2 ' x 6 s l ) + Vl s 2 
+ k 2 s 3 W 3 W l - VsVs + ° 3 T 1 - V 2 

- ZA^Jti^ + ^2 ) C 2 S 2 C 3 S 3^ / C 2 
x 3 = X 4 

x k = X 2 X 6 C 2 ~ V X 2 C 1 + X 6 C 1 S 2' ) “ A 2 S 1 " k 2 C 3 W l W l " k l S 3 W 2 W 3 
+ c 3 v 2 - b 5 y x + 3A 3 (k 2 c| - k x 4 )c 2 S 2 


x 5 x 6 

X 6 = ' X 2 X 4 C 2 ■ X 2 S 2 + A 1^?2 S 1 C 2 + X 4 C 1 S 2^ 

- A 2 C 1 C 2 " Ws + V 3 - 

where = C^/n i = 1*2., 3* 

= e/n = [1 + 2 e cos(t + 0 Q )] 

Ag = 0/n 2 = -2e[sin(* + 9 q ) + 5e cos(t + 0 o )sxh(t + 0 Q )] 


(2.14) 


A = (n/r 3 n 2 ) = [1 + 3e cos(t + 0 O )3 
o 

W x = a^/n = x 2 c 2 c 3 + + A^B^Sg - c^Cj) 

W 2 = co 2 /n = x 4 c 3 - x 2 e 2 s 3 + A^SjCg + c-jSgSj) 
W 3 = a) 3 /n = x 6 + x 2 s 2 + SLffz - 


' (2.15) 


The acquisition equations of motion for satellite (4) are 
obtained by adding 

RHSi = rh S;L /n 2 = A 4 (^ Sl c 2 + 

(continued) 
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(2.16) 


p 

RHSg - rhs^/n = 0 

RHS 3 = rhs 3 /n 2 = \(Vl S 3 - ^c^Kl" 1 ) 

■where = (C D p^Sh 2 /2a 2 ) [1 + 2e cos(t + & Q ) ]exp{Ke[cos(T + 0 q )- 1]} to 
the right-hand side of the second., fourth and sixth of equations (2 13), 
respectively. Equations (2.16) are obtained from equations (2 11a)- 
(2 11c) and (2 12) 

2 The Simplified Steady-State Equation of Motion 

High-accuracy earth-pointing motion or steady-state motion is 
defined to be motion such that |x | < 10 ^ radians, i = 1, ,6 (see 

equations (2 13)) (In the search for a steady-state control law, i e , 
for the functions v = v (x.., • .,x^,t), i = 1,2,3, it is required 

1 |. 3 _ X O 

that |x | S 1 1 X 10~ radians ) In the steady-state mode the terms 
an the nonlinear equations of motion which contain products of some of 
x , i=l,. ,6, are insignificantly small and can be omitted The 

terms m the earth oblateness part of the gravitation torque (terms with 
J as a coefficient) cannot m general be neglected as was done m the 
acquisition equations. However, m the equation corresponding to the 
Tgn 3 component of the torque the oblateness terms can be neglected 
since for an entire orbit these terms are about one one -hundredth of 
the inertia torque term, 0 In the two equations corresponding to the 
T-j n^ and -^—2 com P° nen ^ s °f the torque the oblateness terms which 
are significant are of the same order of magnitude as the largest of the 
other terms for some orbits, but, these oblateness terms which are 
significant for some orbits are insignificant for other orbits (These 
significant oblateness terns are periodic with zero occuring at the time 
of coincidence of B* with a point of the earth's equitonal plane 
except m the case where the orbit is m the earth's equatorial plane 
and all of the oblateness terms are zero.) It should be remarked that 
for greater earth-pointing accuracy, say an order of magnitude greater, 
some of the oblateness terms are the most significant terms m the 
equations corresponding to T^n ^ and for most orbits if the 

aerodynamic torque is insignificant. 
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In the simplified steady-state equations of motion, which are 
written below, the significant oblateness terms are included; -although, 
m some of the control law analysis the orbits are chosen so that these 
terms are zero. Before the equations are written it should be remarked 
that the "best" functions v_^ = v i (x^, .. ,Xg,r), i = 1,2,3, are such 
that v^ are of the same order of magnitude as x , 1 = 1, . ,6, so 
that terms m the equations which are products of v and x are 
insignificant . 

If and c^ are replaced by x and 1, respectively, 

then for satellites (l),(2) and (3) the simplified steady-state equations 
of motion are 

X i = X 2 

*2 - - k l4 X l + A 2 X 3 + W - a iV ( V r)2 V 8 C « + T 1 

X 3 = X 4 

X 4 = “ A 2 X 1 " * k 2 (3A 3 + 4 )X 3 

t 8k 2 A 5 J(r E /r) 2 s 5 c 5 s 0 + v g 

x 5 = X 6 

x 6 = - 3k 3 A 3 X 5 - A 2 + v 3 (2 17) 

where = 1 - k^, = 1 + k^ and A l = 1,2,3, are given m (2.15). 

For satellite (4) the simplified steady-state equations of 
motion are 



x 2 K l A l x 4 + ( a 1)./^i^ 3 + v i 

x 3 = x 4 

(continued) 
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2 

= ^ 2 ^ 1*2 + 8k 2 A 3 J ^ r E //r ' > s 5 c S s 0 + v 2 * 

X 5 = X 6 

*6 - - VA + t 3 (2 l8) 

where A^ is the same as m (2.16) 

E. THE ERROR IN THE MOTION DUE TO THE SIMPLIFICATIONS 

An upper "bound on the motion error due to the simplifications made 
m the equations of motion can be easily obtained for satellites which 
are stable m the sense of DeBra [11]. Suppose that the full nonlinear 
equations of motion are written m matrix form as 

x' = A(t)x + Bv(t) + S 1 [*^v(t),t] + (2.19) 

where the transpose of x is given by x = (x-^x^, . ,Xg), A(t) is a 

six-by-six periodically time-varying matrix,, B is a six-by-three 

constant matrix, v (t) = [v..(t), ,v (t)] is the control vector which 

1 x o 

is a known function of r and the steady-state (or acquisition) solutions 
of (2.19), say cg(T), for given initial conditions, corresponding 

to t , g [x, v(t), t] is a vector function which contains all of the 
forcing function (or nonlinear) terms retained m the simplified steady- 
state (or acquisition) equations, and, gg[x,v(T)] is a vector function 
which contains all of the terms omitted m the simplification Then the 
solution of (2.19) can be written as 



T 

+ J ®(T,A)j^[cp(A),v(A),A]aA (2 20) 

T 

o 
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where is the fundamental matrix of x* = A(x)x. 

If the norm is defined by 

6 6 

Ixi = ^ |xj and |$j |«P XJ I 

1=1 i,j=l 

where 9 are elements of $ and if 9 (t) is defined to be the 
1 J s 

solution of the simplified equations, then 

T 

|<9 (t) - ^(t)! 5 J |$(t,A) | [g g [9 (A),v(A),A] |dA (2 2l) 

T 

O 

m the steady-state case since g^ is a function of t only m this 
case 

DeBra has investigated the motion of satellites m elliptic orbits 
about an oblate earth and has found for certain satellite configurations 
that the motion is "stable" Satellites which have such "stable" motion 
have the characteristic that m ? X A) j & 20 for t ^ A ^ t In 

the steady-state case, since v(t) is such that |9 (x) | si 1 X 10~ ' 
for all t, the maximum value of I 1S ^ ess than 10 

Since the equations are periodic with period 2at, the equations for the 
steady-state case need be integrated only over the interval x„ - x = 2it. 

/* X o 

Thus, m this case |9 (t„) - 9 (x„) | < 1 26 X 10” Since j9 (x„) | « 

— x X si 1 

10 ,i=l, ,6, for most solutions obtained (see Chapters IV and VI ), 

the error in 9 (t) is generally less than 2 $ for satellites which are 

S X 

stable m the sense of DeBra Satellites (l) and (4), the roll and yaw 
motions of satellite (2) and the roll and pitch motions of satellite (3) 
are "stable" 

Since the vector function g^ is a function of cp(x) as well as 
of x m the acquisition case, no meaningful upper bound on the error 
due to simplifying the acquisition equations can be found with the above 
method However, Busch [6] found that there was little detectable error 
m his acquisition motion obtained from his simplified equations, which 
were completely linearized. 

Ji Busch compared the linear acquisition solutions with the nonlinear 
solutions 
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F. DISTURBANCES 


Those influences of the motion (or those torques) which are not 

taken account of by the equations of motion axe called disturbances The 

primary sources of possible disturbances are the sun,, the earth's magnetic 

field and meteoroids. (Other disturbances, eg, misalignment of gas 

jets, are discussed m Chapter VI ) 

The torque for due to the sun's gravitational attraction of the 

-4 

satellite is no larger than 10 times the torque, T , due to the 

& 

gravitational attraction of the earth when the satellite is m the high- 

accuracy earth -point mg mode of motion This torque is, of course, very 

insignificant when studying the motion for one orbit or less. 

Emissions from the sun exert a pressure on the exposed part of the 

satellite's surface area. A torque for B* due to this pressure (usually 

called solar radiation pressure) can be quite significant Expressions 

for this torque have been derived by McElvam [26] and Wheeler [34], 

From these expressions it is concluded for the satellites considered here 

that the solar pressure torque is no larger than 10” times |T j when 

the satellite has large motion and a distance, d , between B* and the 

s 

center of solar pressure of 1 ft For high-accuracy earth -pointing 

motion the solar pressure torque can be as large as |t | if d is of 

g s 

the order of 1 ft This torque can affect the form of the required 
steady-state control law and the cost of the control. However, since 
the solar torque is very similar to the aerodynamic torque m effect and 
since the aerodynamic torque is accounted for, the solar pressure torque 
per se is not considered further 

The torque for B* due to the interaction of a satellite with the 
earth's magnetic field can have a significant effect on the motion 
Bandeen and Manger [2] have considered a model of the Tiros I satellite 
for correlating data on the precession of the nearly earth-pointing spin 
axis. The model included the effect of the interaction with the earth's 
magnetic field To model this satellite Bandeen and Manger assumed that 
a circular conductor with a one meter diameter was on board and carrying 
a current of one ampere. This model correlated well with the data 
received. The magnetic torque of this model was of the same order of 
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magnitude as the earth's gravitational torque for an inaccuracy m 
earth-pointing of 0 1 radians However., since Tiros I was much less 
massive and extensive than the satellites considered here (with about 
the same asymmetry) , the same magnetic torque is only about 10 times 

|t | for the same earth -pointing accuracy. If this same current- 

— g 

carrying coil were placed on the satellites considered here., the magnetic . 

torque would be of the same order of magnitude as T when the satellites 

— g 

are m the high-accuracy earth-pointing mode of motion. 

For actual satellites m general it is difficult to say what effect 
the torque due to the interaction of a satellite with the earth’s 
magnetic field will have on the motion Part of this torque results from 
a residual component due to magnetization of some parts of the satellite 
Since the magnetic torque is considered to be no larger than the other 
torques and since by properly designing the satellite this torque can be 
made quite small, it is considered hereafter only as an unknown distur- 
bance of a certain maximum amplitude which must be overcome by the « 
control 

The forces and their torques due to collisions with meteoroids are 
considered to be extremely insignificant for nearly the entire lifetime 
of the satellite Christman and McMillan [8] have concluded from 
measurements obtained by the Explorer, Mariner and Pegasus satellites 

that the probability of no impact by a meteoroid with momentum of 

-3 

6.67 x 10 lb-ft/sec is 0 99 and that for meteoroids with greater 

momentum impact it is even less. Using experimentally available data, 

Whipple [35] has estimated that meteoroids with speeds of 20,000 ft/sec 

-12 

and masses of 2 x 10 lb. will imp act. a satellite with a cross- 
sectional area of twenty square feet about every 20 seconds at an 
altitude of 200 miles. Whipple has also concluded that the frequency 
of impact decreases log arithmetically with altitude and with increase 
m the mass of the meteoroid 

Cloutier [9] has found that the torque due to the forces exerted on 
a satellite by meteoroid impact can be as great or greater than the 
torque due to the earth's gravitational attraction but that this Is 
generally a rarity A provision must be made for such rare occurancesj 
and, this is done m Chapter IV. 
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A summary of the most significant torques, which can he predicted 
with some certainty, is found m Table (2 2) In this table for various 
ranges of 1 = 1*2,3, the satellite numbers ((l) thru (4)) are 

placed under the torques and beside the component designation if for a 
particular satellite the component of the torque is significant A 
component of a torque is insignificant if it is at least one order of 
magnitude less than corresponding components of other torques For 
satellite (l) some components of the totality of the torques m the 
table are shown as insignificant for some ranges of |0 |, 1 = 1*2,3. 

The reason for this is that these components are insignificant compared 
to the corresponding components of the inertia torque, which is not 
included m the table. 

The model used for determining the magnetic torque m Table 2.2 was 
the same as Bandeen and Manger used for Tiros I The expression used 
for the solar pressure torque of Table 2 2 was derived by Wheeler. 
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TABLE 2 2 THE SIGNIFICANT TORQUE TERMS FOR VARIOUS EARTH-POINTING ACCURACIES 
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Ill THEORY OF THE CONTROL OF THE SATELLITES' ATTITUDES 


A THE MOTION REQUIREMENTS IN TERMS OF THE STATE OF THE PLANT 

In Chapter 1, Section A, the attitude motion requirements of the 
satellites -were given m geometric terms Here., they will be given m 
an algebraic form 

In Chapter II the differential equations of attitude motion were 
derived m terms of the three -axes Euler angle s, which define 

the orientation of the satellite with respect to the orbiting earth- 
pointing reference frame On page lij- of Chapter II the three second-order 
differential equations of attitude motion were replaced with six first- 

*b 

order differential equations by defining the vector x = (x^,Xg, ,x^) 

The state of the plant at some time,, t, is defined to be the six- 
dimensional vector x(t) = [x^t^x^t), ^^(t)]^ The state space is 

the six-dimensional euclidian vector space, X, of which x(-r) is an 
element 

High accuracy earth -pointing requires that x , i = 1,2, ,6, be 

kept less than or equal to given small positive numbers for the lifetime 
of the satellite. (See Appendix D ) Let s , l - 1,2, 6, be these 

given* small positive numbers Then it is required that x^ ^ s^, 
l = 1,2, ,6 These inequalities define a closed and bounded region, S, 

of the vector space X which contains the origin or zero vector 

B SATISFACTORY PERFORMANCE AND THE COST 

A satellite attitude controller (see Figure 3 l) will he said to have 
satisfactory performance if acquisition to the region S from x pa 1 50 
< tc/2, i = 1,2, ,6, can be accomplished with near minimum cost m less 

than the time of one quarter orbit and if the station-keeping part of the 
controller with near minimum cost for the lifetime of the satellite can 
keep the state space trajectory from departing the region S by a signi- 
ficant amount (except when large unaccounted for disturbances overpower 
the station-keeping part of the controller) 

~ 

See Chapter IV 
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Figure 31 A Simplified Diagram 












The cost of a satisfactorily performing satellite is ultimately a 
combination of the monetary cost and time costs of producing the satellite 
which performs the required mission. In this presentation it is assumed 
that except for the attitude controller the satellite is of a given con- 
figuration and cost., so that, reductions m cost are obtainable only by 
reducing the cost of the controller 

The cost of the controller is not just the cost of constructing the 
controller and operating it m orbit but includes the cost of placing 
the controller m orbit with the satellite. A heavy controller will 
require more power from the vehicle which orbits the satellite than a 
light controller 

When the controller uses gas jets to provide the control torque, the 
weight of the fuel for the gas jets can increase the cost even if increases 
m the cost of orbiting the satellite are not considered For example, 
the increase m cost due to an increase m the weight of the fuel can 
result from expensive packaging caused by the increase m volume of the 
fuel container or from the increase m the cost of constructing the fuel 
container which must withstand higher pressures if the volume is kept 
small 

The feedback attitude controller will contain electronic computing 
elements and gas jet thrusters The number, sizes, weights and complex- 
ities of these should he kept at a minimum for minimum overall cost It 
is true that the use of off-the-shelf hardware components will reduce the 
cost if the weights and complexities of these components are not pro- 
hibitive . 

Thus, to reduce the cost of the attitude controller the weight of 
the controller should he reduced as much as possible compatible with 
inexpensive off-the-shelf components and with the simplicity of the over- 
all system 

The power supply for a year or more of control is the heaviest com- 
ponent of the controller so that the weight of the power supply should he 
at the focus of attention. Other components are nearly fixed m weight 
by the state of the art except for the thrusters which generally decrease 
m size and weight with decrease m trust magnitudes It should he 
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pointed out, however, that for very small thrust requirements hot gas 
or ion thrusters are used and they require a weightly power supply for 
supplying heat 

For a controller which utilizes gas jets for control the power 
supply consists of the fuel, which is usually an inert gas, the fuel 
container, tubes for carrying the fuel to the thrusters and pressure 
regulating devices Except for the fuel, the components of the power 
supply are nearly standard m weight for all but high gas pressures By 
reducing the weight of the fuel, the gas pressure and/or volume will be 
reduced Thus, the weight of the fuel will be at the focus of further 
attention 

C THE OBJECTIVE 

In summary, the objective is to determine a controller which will 
acquire the region S of the state space within the one quarter orbit 
time limit from values of x , i = 1,2, ,6, of about 1 5 radians and 

which will keep the state of the system within S for the remaining 
lifetime of the satellite while using as little fuel as possible compatible 
with costs due to the complexity of components and the development of new 
components 

D. MATHEMATICAL TOOLS USED FOR DETERMINING A CONTROL LAW 

Two basic approaches to solving the problem of determining a control 
law are 

1) The application of mathematical results based on fundamental 
principles which apply to the minimization problem 

2) The simulation on an analog or digital computer of the plant 
with various controllers which are determined by theoretical 
knowledge of the behavior of the plant under the action of the 
controllers 

In the mathematical theory the problem to be solved is one m the 
field of optimal control with inequality constraints on the state variables 
(the station-keeping part) and without such constraints (the acquisition 
part) 
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Pontryagm's Maximum Principle [29] and extensions of it will be 
used as a mathematical aid m obtaining solutions to both the acquisition 
problem and the station-keeping problem. For future use Pontryagm's 
Maximum Principle is stated below m the form which is applicable to both 
problems . 


1. Pontryagm's Maximum Principle for Ufonautonomous Systems 

The maximum principle gives necessary conditions for the solution 
of the optimal problem 

Suppose the following are given 

1) The differential equations of motion 

x T = f(x,v, t) (3 l) 

where x = (x^ ,x n )^ and f ~ ^ • jf 11 )^* (For example, 

equation (3.l) can represent equations (2 It).) 

2) An initial point m the state space, say x^, which describes 

the motion at an initial time, r , after which the motion is 

o 

considered to he under the influence of the control, v 

3) A final point m the state space, say x^,, which describes 
the motion at the final instant of consideration, 

k) The cost functional 


J = J f°[x(r), v( t ) ] dr 

T 

O 


(3.2) 


Then the optimal problem is to find the control v = v(t), which 
is at least as smooth as piecewise continuous, that causes the motion to 
be given by x f at and results m the minimization of J 

It should he pointed out that the maximum principle applies to 
problems in which the control is much less smooth than piecewise continuous 
However, there is no need to consider a more general class of control since 
the control devices considered here are adequately described by piecewise 
continuous functions 

The set V is defined as the set of all bounded piecewise 
continuous vector functions v = v(t) The function H[p(r), x(t), t. 
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v(t)], the Hamiltonian; is defined by 

n 

H(P;X;T;V) = ^ P^Cx^VjT) (3 3) 

1=0 

"t 

where p = (p A ;P-,; .,p ) is defined m detail below H . (p;X;ir) is 

u j. xi max ~~ 

given by 

H v [p(t), x(t); t] = max H[p(-r); x(t); t, v(t)] (3 4) 

maX “ v( t ) eV 

where the maximization of H on v(t) is with respect to p(t); x(t) 
and t 

In order that v(t) yield a solution of the given optimal 
problem it is necessary that there exist a nonzero continuous vector 
function p(t) = [pq(t); p^(t); •;P n ('t')]^ corresponding to the functions 

v(t) and x(t) through equation (3 l) and 

P ' = - cHl/S X^; 1 = 1 ; ,EL (35) 

such that 

( 1 ) for all T; ^ t s t , the function H[p(x); x(t); t, v] 
of the variable ve? attains its maximum value at v = v(t), 1 e ; 

H[p(t), x(t); t, v(t)3 = H Cp(t); x(t); t] 

lifCLA. *” 

where H is given by equation (3 3) and H is given by equation (34); 

max 

( 11 ) the function p (t) 1S a nonpositive constant 

This is a statement of Pontryagm's Maximum Principle (PMP) for 
nonautonomous systems. Pontryagin, et al [29] assume that f°(x; v) 
and f(x; V; t) have continuous first derivatives m x in their proof 
of the maximum principle This restriction m the proof can be weakened 
so that the maximum principle applies to problems m which 
v(t))/ 3 x is only continuous almost everywhere. Breakwell [3] gives a 
derivation m this case for autonomous systems m an early paper A slight 
extension of Halkm's work [17] results m a proof Eozoroer [30] gives 
a proof for optimal problems m whach both f°(x; v) and f (x; V; t) have 
continuous second derivatives m x This proof can be modified so that 
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it is clear that PMP applies to problems which have f°(x, v) smooth 
only almost everywhere along x(r), s t £ 

The maximum principle as stated above can be extended to apply 
to the optimal station-keeping problem and to the optimal acquisition 
problem m the case acquisition to the region S is required rather than 
the case were acquisition to the origin of the state space., X, is 
required The applications of the maximum principle to these problems 
differ primarily m the boundary conditions 

Before the maximum principle is applied to an optimal problem 
there should be some certainty that the results of the application will 
give correct information about the solution of the problem The maximum 
principle gives conditions which, if they are not satisfied, imply that 
the control is not optimal. 

Suppose that there are many controls veV which take x to 

x f and suppose that only one of these satisfies the maximum principle 

Then this one control is the optimal control if an optimal control exists 

An optimal solution does not exist if the functional J for the control 

veV which takes x to x^ while satisfying the maximum principle is 
“ — o — x 

not a minimum. Otherwise, a solution to the optimal problem exists. 

In the present investigation Pontryagm’s Maximum Principle 
gives a complete set of relations for the determination of a control 
Even so and even if a solution exists, there is no guarantee that the 
control obtained from the maximum principle does not give a local minimum 
of J 

In the linearized acquisition problem an optimal solution exists 
for veV and the functional J does not have local minima, so that for 
a solution of the optimal linearized acquisition problem, Pontragm's 
Maximum Principle gives sufficiency conditions. A proof of this can be 
found m Rozonoer [30] for the case when f(x,v, v) and f'°(x, v) have 
continuous first derivatives m x which is true m the acquisition 
problem 

In the application of the maximum principle to the approximate 
station-keeping problems, the function f°(x,v) is either nonlinear or 
has derivatives m x which are only almost everywhere continuous m 
x(t). In this case Rozonoer* s proof of the sufficiency of the conditions 
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of the maximum principle does not hold. However, since the conditions 
of the maximum principle are necessary conditions and since they are 
sufficient to determine a unique motion, then they are sufficiency con- 
ditions for the optimal solution if an optimal solution exists In the 
case of the approximate station-keeping problems, the existence is 
concluded by direct reasoning 


2. Application to the Acquisition Problem 


Busch [6] has found a nearly optimal feedback control which will 
cause the motion to proceed from x to near zero for a "stable" satellite 
configuration For simplicity and reliability Busch's control law is a 
function of the state only The method Busch uses for determining a 
control law from the maximum principle is reverse-time integration, i.e , 
once v = x(p;Xj> T ) ■* ias ” been found from the maximization of the Hamiltonian 
H (see equation (3 4)), the equations (3 l) and (3.5) are integrated 
backwards from p(tu) and x„ to x and p(r ) The solutions x(t) 
and p(t), t q s t ^ are then analyzed m order that characteristics 
of the solutions x(t) and p(t) can be found which enable the control 
to be written as 


v = v(x,t) (3 6) 

In the reverse time integration procedure T is replaced with t* = T f “ T 
and x is determined by the choice of p,, and the interval of mtegra- 
tion, - t (when x f is given) 

In Chapter V a solution to the problem of the optimal acquisition 
control of unstable satellites is given This solution is based on 
Busch's solution, the maximum principle and the imposed time limit of 
acquisition 


3 Application to Approximate Station-Keeping Problems 

Approximate solutions of the station-keeping problem can be 
obtained from PMP by taking the cost functional to be of the form 



T 

o 


(3 7) 
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where g(v) is scalar function of the control vector and f (x); a 
"penalty function " , is a scalar function of the state vector. 

Since the minimization of the weight of the fuel is required; 
g(v) is chosen as 
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i=l 


( 3 . 8 ) 


The optimal control for the problem with the cost functional given by 
(37) with (3 8) and any f(x) corresponds to the use of simple gas jets 
The "penalty function" is chosen such that the control keeps x m 
or very near S while using the least possible amount of fuel. If 
f(x) is chosen so that it is zero when x is m S; the functional 
(3.8) to be minimized becomes the minimum fuel functional while x is m 
S Thus ; it seems that the fuel expenditure should be a minimum at least 

for those period when x is m S 

Possible nonnegative functions which are zero when x is in S 

are 


r 


if x ^ s 
1 1 1 1 


6 

Y — 1 ^ 

(5) = ) < 10 1 ( x 1 - s^) if \ > S 1 


n 


-10 1 (x + s ) if x < - s 


1 1 


( 3 . 9 ) 




i=l 


1 
k. 2 


n 


if x ^ s 
1 1 ' 1 


10 1 (x -s ) if x > s 
11 11 


n 


10 1 (x +s ) if x < - s 
1 1 1 1 


(3.10) 


(continued) 
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r 


6 n 

f g (x) = Y < exp [ 10 1 ( V b i )]-1 if x^ > s x 


if X \ £ s 

i X 1 ! 


(3 ll) 


L- 

1=1 


n 




exp [-10 1 (x i +s i )]~l if x i < - s^ 


where n , i = 1, . ,6, are numbers chosen so that x stays m or near 


Functions similar to 


6 


,w -y \ 

1=1 


Ki/<h 


s 


Kb 


(3 12) 


where k^ and (a number slightly greater than one) are chosen to 

keep x near S, are possible choices Since they exhibit singular 
behavior as |x^ | -* & • s , they result m the state staying very near 
the boundary of S or in S (depending on £ ), but, since they are not 
zero for x m S, fuel is wasted 

For an idea of the relative values of the functions given by 
equations (3 9), (3 10), (3.1l) and (3 12) for values of x see 
Figure 3.2. 

Since m the station-keeping problem x^ and x f are arbitrary 
to within being m S or on the boundary of S, the boundary conditions 
on the variable p(r) are somewhat better known m advance m the station- 
keeping problem than m the acquisition problem 

In Pontryagm, et al [29] and m Eozonoer [30] it is shown that 
a necessary condition for an optimal solution is 

£(T f ) = Cp 1 (' r f )^ • ^P n ( T f )] t = (°> #°) t ( 3 13 ) 

if x(T f ) is (free) m the interior of S or 

p(T f ) = [p 1 (T f ), ^P n (T f )3 t = - ll b(x f ) (3 Ik) 

if x(t f ) is (free) on the boundary of S She vector b(x f ) is the 

outer normal to the boundary of S at x( ) , p is a nonnegative constant 
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and n -6 in the present investigation 

The relation expressed by equation (3 Ik) is called the trans- 
versality condition Since S is defined by |x | £ s , i = 1, ,6, 

equation (3 lk) can be written as 

P 1 (T f ) = - il sgn [x i (T f )], 1 = 1, ,6 (3 lk 1 ) 

The maximum principle requires that p Q be a nonpositive 
constant Rozonoer has shown that if p = - A, then a necessary con- 
dition is A > 0 Since H is homogeneous m p , i = 1, ,6, the 

maximinization of H is independent of A For convenience choose 

A = 1 

The initial value P.( T 0 ) 1S m general not known If x(t q ) is 
required to be on the boundary of S, then ^(t ) must satisfy a con- 
dition similar to that given m equation (3 lk) Otherwise, p(t q ) ls 
unknown Since x(t ) is not required to be on the boundary of S, the 
transversality condition for x(t ) is not applied 

It should be noted that the boundary conditions given by equation 
(3 lk) do not apply when x is on a "corner" of S, i e , a point on 
the boundary of S at which two of ]x.J, i = 1, ,k, or both of \x. | 

and jxg [ are equal to the corresponding numbers s This follows from 
the fact that the normal to the boundary of S is not defined at a 
"corner" If this should present a problem, the "corner" can be smoothed 
out by constructing a suitably smoothe "surface" m the "corner" between 
the faces of the polyhedron m the six dimensional space (if these 
surfaces are small enough, no physically measurable changes m the boundary 
will occur ) 

The above boundary conditions and the cost functionals given by 
(3 7)-(3 11) are used m Chapter IV with the maximum principle to deter- 
mine approximate solutions of the steady-state problem 

Similar approximate methods of solution for linear, minimum time 
problems with restricted phase coordinates has been developed by Lee [2k] 
and by Russell [31] Lee applies the method to the minimum time acquisi- 
tion problem of a l/s plant with bounded phase coordinates 
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4 Necessary Conditions for an Exact Solution ( I'JCES ) of the 
Station-Keeping Problem 

Several investigators have to some extent developed exact solu- 
tions to the hounded phase coordinate optimal control problem Bryson, 
Denham and Dreyfus [5] have found necessary conditions for solving 
problems m which the control is a scalar and certain smoothness assump- 
tions must hold Chang [7] and Russell [31] are concerned with sufficiency 
conditions for solving the linear, minimum time acquisition problem with 
bounded phase coordinates. Pontryagm, et al [29] devote their Chapter 
VX to obtaining necessary conditions for the solution of the optimal 
control problem with restricted phase coordinates These necessary 
conditions, which are discussed m this section, will be used m the 
derivation of the station-keeping control law m Chapter IV 

Since the control function v(r) is piecewise smooth (discon- 
tinuities of the first kind) and since the other terms m the right hand 
side of the differential equations of motion are smooth, there exists 
only a finite number of points m time at which the trajectory, x(t), 
is on the boundary and either x(t ) or x(t') or both are not on the 
boundary If x(t) 1S ' tlie state at such a point m time t and if both 
x(t ) and x(t + ) are m the interior of S or one is on the boundary 
of S, then the point x(r) is called a junction point of the trajectory 
Of course, there are only a finite number of junction points 

In the derivation of the necessary conditions for an optimal 
trajectory with restricted states it is assumed that* (l) the optimal 
trajectory has only a finite number of junction points, (2) the optimal 
trajectory lies either entirely m S or on the boundary of S, and, 

(3) the parts of the optimal trajectory which lie on the boundary of S 
for a finite interval of time must be "regular" 

The boundary of the region S must be somewhat "regular" or 
smooth, 1 e , if s(x) = 0 describes the boundary of S, then grad s(x) 
must be continuous and not vanish for any x m the boundary of 
S In the station-keeping problem with S defined by jx^J s s , 

1 = 1,2, ,6, this is not the case unless the "corners" of S are 

avoided when the trajectory coincides with the boundary of S By 
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suitably "smoothing" the "corners" (see p 35 of this section) this 
mathematical difficulty is overcome In Chapter IV the HCES are applied 
to several approximate single -axis satellite motions by replacing the 
"non-regular" boundaries of a phase-plane projection of S "with appro- 
priate "regular" boundaries These "regular" boundaries are section of 
known optimal trajectories which must exist and deviate from the phase- 
plane projection of S (for bounded v) by some small allowable amount 
In Chapter VI of Pontryagm, et al ; the concept of regularity 
is used m the derivation of necessary conditions for optimality of those 
parts of the trajectory (if any) which lie entirely on the boundary of S 
It should be noted that Pontryagm, et al , derived these necessary con- 
ditions under the assumption that each of f°(x,v) and f (x, v, t) have 
continuous first derivatives m both x and v If it is assumed that 
f°(x,v) has only continuous first derivatives m x and not m v, then 
a weaker set of necessary conditions are obtained These weaker conditions 
are the same as Theorem 22 of Pontryagm, et al , except that dH/dy 
does not exist everywere along the trajectory which coincides with the 
boundary 

If there exists an optimal trajectory which lies m the region 
S, it is possible that part of this trajector lies entirely m the 
interior of S These parts of the optimal trajectory must satisfy the 
conditions of the maximum principle 

If the regular optimal trajectory contains only a finite number 
of junction points, then an additional condition on the function p(t) 
at the junction time can be derived Pontryagm, et al , call this 
condition the jump condition The jump condition is satisfied if one of 
the two following conditions is satisfied 

£(t + ) = p(t“) + p grad s[x(t)] 

E(t + ) = p(t”) + P grad s[x(t)] = 0, p f 0 (3 15) 

where p is a number to be determined 

(in the station-keeping problem the region S has corners In 
a note m Chapter VI of Pontryagm, et al , (page 310) it is pointed out 
that jump conditions completely analogous to those of equations (3 15 ) 
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must be satisfied at the transition point from one smooth part of the 
boundary of S to another ) 

The necessary conditions for an exact solution of the station- 
keeping problem can be summarized as follows- The conditions of the 
maximum principle must he satisfied by each part of the regular optimal 
trajectory which lies m the interior of the region S The conditions 
of Theorem 22 of Pont ry agin, et al , must he satisfied by those parts of 
the optimal trajectory which lie entirely on the boundary of S The 
jump condition must be satisfied at a junction point 

The above necessary conditions are, generally, insufficient to 
determine the optimal control Without any conditions other than those 
above the search for the optimal control must, usually, he carried out 
m a large dimension parameter space Parameters of this space include 
the number of junction points, the number of the parts of the trajectory 
which lie on the boundary, the boundary conditions p( T 0 ) and p(^) 
and the number (i 

In the next chapter the above necessary conditions will he used 
m a search for a control law for the station-keeping problem This 
search will also be aided by the optimal solutions to the several approxi- 
mate station-keeping problems 
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IV DERIVATION OF STATION-KEEPING CONTROLS LAWS 


The attitude motion of a satellite m the steady-state mode of 
controlled motion is to be such that ♦ 

1) The initial instant of station-keeping, t o , is arbitrary, 

1 e , t q is the time at which the satellite is at an arbitrary 
point in its orbit, 

2) The motion at time t q is arbitrary to within x(t ) being m 
the region S of the state space, 

3) The motion after r and until some given final time t_ must 

o f 

be the result of a control v(t), s -r ^ t^, which keeps 
x(t), t o < t £ r f from departing the region S and which uses 
as little fuel as possible 

In the application of the theory of the previous chapter to the 
search for station-keeping control laws, it is convenient to take = 0, 
T f = 2rt and (m either (2 17 ) or (2 l8)) Q q = 0 This can be done with- 
out loss of generality since (l) x(t q ) and x(t„) are arbitrary to 
within their being m S (x(t ) can be the final state of an acquisition 
trajectory or the final state of a previously considered station-keeping 
trajectory), ( 2 ) The equations of motion are periodic with a period of 
2it radians, and, (3) In the lifetime of a typical satellite the boundary 
of S is encountered many thousands of times and the region S is 
nearly covered by the state- space trajectory (see Part 1 of Section A) 

The region S was defined to within the numbers s , 1 = 1, 6, 

m Section A of Chapter III If s , 1 = 1, . ,6, are greater than 

-2 1 

about 10 , the simplified steady-state equations of motion do not give a 

suitable description of the motion The lower limits on s , 1 = 1, ,6, 

are fixed by the state-of-the-art m the construction of sensors and 

controllers which have very small gas jet thrust and/or time delays 

- -if. 

Hereafter, the numbers s , 1 = 1, ,6, are assumed to be 10 unless 

otherwise specified 
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A THE APPROXIMATE SOLUTIONS - SPECIAL ORBITS 


In this section several approximate solutions will be presented 
They are approximate solutions to the problem of determining a control 
which will keep the state of a satellite 1 s attitude motion m the region 
S while minimizing the fuel used 

The approximate solutions are approximations m the sense that the 
state space trajectories are allowed to exist the region S by a small 
distance (compared to the maximum dimension of S) and m the sense that 
integral constraints on the states are used to limit the motion If the 
integral constraint is given by 

T f 

Os J f[x(t)]dT s A 

T 

0 

where A is a given (perhaps small) positive number and f (x) is given 
by (3*9) or (3 10) or (3 ll), then this optimal problem is equivalent 
to the optimal problem of Part 3, Section D of Chapter III The problems 
are equivalent m the sense that their solutions (as obtained from the 
maximum principle) are the same 

In this section the equations of motion are assumed to be given by 

(2.17) with the angle S either zero or jt/2 radians These values 

correspond to satellites m nearly equatorial or nearly polar orbits, 

respectively They are used initially (Part l) to simplify the analysis 

although the results do not depend on 5 In Parts 2-b these orbits, 

i.e , these values of & are assumed so that equilibrium points exist 

In Parts 2, 3 and 4- of this section the equations of motion are 

approximated further The "stable" single-axis motions of the satellites 

2 

are approximated by x + a x = v, and, the "unstable" single-axis motions 

2 2 

are approximated by x" = v and x" + a x = v with a <0 The study 
of these simple motions results m characteristics of a suboptimal 
minimum-fuel station-keeping control as well as a check on the methods 
of Part 1 
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1 A Station-Keeping Control Law Obtained From PMP with 


rf 

J =J [^(x) + 

T 

O 


I 


|v i I ]dT 


The Hamiltonian for the system of equations (2 17) with 6=0 
or 6 = tc/ 2 can be written as 

H = P 2 v x + p^ 2 4 pgv 3 - |v 1 | - |v 2 l - |v 3 | + (terms which do not 
contain -v > l = 1,2,3) (4 l) 


since as was seen m Chapter III p^ can be taken as -1 

The control v cannot be optimal unless it maximizes H for 
all values of p, x and t Thus, from (4 l) the optimal control for 
bounded v must be the "coast function" of p given by 


v 

i 


= CST(p 


2 x 1 ) 



sg H (p 2ii ) 


If 

If 


IP. 


2X1 


2x1 


< 1 
£ 1 


i 


1 , 2, 3 


where N , i = 1,2,3, are given positive numbers and SGW(g) = (g)/|(g)| 
Since the last two equations of (2 17), the pitch equations, 
are not coupled to the first four of (2 17 ), "the yaw -roll equations, the 
controlled pitch motion can be solved for independently of the yaw-roll 
motion These equations (for the cases when 6 « 0 or 6 it/2 ) 
with the corresponding equations for the adjoint variables as determined 
from (3 5) can be written m backwards times by letting t* = 2it - t as 
follows 


yaw -roll 


x, = - x„ 


'*2 ‘ Vl*l + *2*3 ' K 1 A 1 3 % - 


x 3 " - K k 


'x^ - - A 2 x 1 + ^2^1 X 2 “ k 2^ 3A 3 + A 1^ x 3 

4l 


- v. 


(4 3) 



'«! ■ - r(x i ) - k i4 p g + Vit 
'Pg * F(x a ) + Pi - K gVit 

'p 3 - - f(x 3 ) ' Va + k g (3A 3 + 

'p^ = F(x 4 ) + K^pg + p 3 (k.k) 

pitch 

' x 5 = " x 6 

' x 6 = ^3^5 " A 2 ' v 3 (4 5) 

' P 5 = - F (x^) - 3k 3 A 3 pg 

*P 6 = F(xg) + P 5 ( k 6) 


■where * () = d( )/dT*, the functions A^, and 
in (2 15 ) hut with t replaced by t* and 



r 0 , 

If 

|x | £ iCf 4 




1 i' 

F(x i )=5f 1 (x)/5x ;L = <j 

io ni , 

If 

-4 

x^ > 10 * , 



If 

-4 

x < - 10 H 
1 


A 3 are the same as 


1 — lj (^7) 


Equations (if- 3)-(4.6) with (4 2 ) and the boundary conditions 
given by (3 13),(3.l4') and x(t*) = x(t^) are sufficient for determining 
the controlled motion, x(t*), the adjoint variables p(-r*) and, hence 
the control v(t*) for -r* £ x* <: t* These solutions can be used to 
determine the optimal feedback control law for this problem or at least 
characteristics of the optimal feedback control law which can be used m 
the construction of a minimum-fuel suboptimal feedback control law 

Although the differential equations are piecewise linear, it is 
not practical to invest a great deal of time m a search for their exact 
solution since the coefficients are time-varying To determine the feed- 
back control law it is only necessary to determine the values of x(t*) 
and T* at which the adjoint- variables assume the values of +1 or -1 
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Thus , the solution of the equations with the aid of a high-speed computer 
is perhaps the most logical approach to determining the feedback control 
law 

After a thorough investigation the hope of using the analog 

computer for solving the above equations was abandoned The reason for 

this was that some of the computer variables., 1 e , the scaled dependent 

4 

variables , are generally 10 times the other computer variables., so that, 
it was impossible to obtain an accurate and therefore meaningful solution 
with this computer 

The equations can be solved accuarately enough with a digital 
computer A Kutta-Merson integration routine was programmed as an ALGOL 
procedure on the Burroughs B5500 computer The Kutta-Merson procedure 
used was a modified version of the Stanford Computation Center Kutta- 
Merson procedure These modifications, the additions of an absolute error 
bound and a step size -cutting limiter, were made to the procedure to reduce 
the computation time The modifications did not affect the accuracy m 
less than the fifth significant figure m the test runs made A further 
reduction m the computation time is accomplished by scaling the equations 
so that the dependent variables are more nearly the same size A listing 
of the program used for integrating the equations is given m Appendix F 

Several computer solutions of the yaw-roll equations and of the 
pitch equations were needed to determine the "best" values of n^, 

1=1, ,6, m the penalty function and H , l = 1,2,3 (These "best" 

values are the values which result m the fuel cost being as small as 
possible while the control keeps |x | ^ 1 1 X 10 \ l = 1, ,6 ) 

The initial choice of values for n and N was made m the 

l l 

following way The possible choices on N lie in a range from the 
smallest values with which control will be maintained at all times up to 
the largest values which cause changes in the motion to occur too rapidly 

In practice the smaller values are the logical ones to choose since with 

large H the lhherient imperfections m the controller can cause unsat- 
isfactory motion and wasted fuel The smallest possible values of N 
such that control could always be maintained, even if |x |, i = 1, ,6, 

grew to values as large as 2 0 x 10”^, where chosen For satellite (2), 
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-4 -4 

for example, these values are ^ = 2 1 X 10 , Ng = 7 1 X 10 and 

U 3 = 2 06 X 10' 2 

The possible values of n_^ are the real numbers From the cost 
functional for this problem it can be seen that its minimum for intervals 
of time when x is not m S is the result of minimizing the part which 
is a functional of x if 

4 2 

^ ± 10 3 -(x i 1 10- 4 ) 

1=1 1 

and 

6 n 

± 10 1 (x i 1 10" 4 ) > N 3 

i=5 

If it is desired to keep |x | ^ 1 1 X 10 , then possible "best" values 

of n^ for satellite ( 2 ) are n^ = 1 0, n^ = 2 0, =10, n^ = 2 0, 

n^ = 2 0 and n^ = 4 0 In the choice of these numbers nearly twice 
the weight was placed on the parts for x^, x^ and since these 

variables change more rapidly than x n , x„ and x R 

The results of a computer solution of the equations for satellite 
( 2 ) with the above values of n^ and are given m graphical form 

m Figure 4 1 From this solution it is obvious that the above values of 
n^ and IT are not the "best" values Even though this solution 
exhibits characteristics which aid both the search for the "best" values 
of n and N and the search for the feedback control law, the cost 
(in computer time) of continuing this method without other aids is 
prohibitive 

An approximate solution of the adjoint equations of this problem 
for arbitrary x(t 75 ') can be obtained so that an optimal feedback control 
law for tms problem can be written m terms of n_^ This approximate 
solution can be used as an aid m the search for the "best" values of 
n and IT and as a guideline m the search for nearly optimal and 
practical feedback controllers of the attitude motions of satellites 
The approximate solution, which is derived m Appendix E, is 
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a Taw Trajectory for = 2 1 x 10 - ^ 



0 0 1 571 3142 4 711 6 282 


TAUSTAR 

b. Taw Adjoint Variables for = 1 = 2 0 

Figure A 1 Solution of Equations (4 3)~(A 6) for Satellite ( 2 ) 
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P x (t*) = P-jgO -[(Bg - Bj^) + (<n 2 + <d 2 ) ](At*) 2 ) 

- P 2e [B l A ' r " S ' + B i i . B 5(^ T+ ) 2 ^ “ P 3e B 5 (AT"0 2 - P^I^At* 

- B^CAt*) 2 ] - P(x 1 )At* - F(x 2 )B 1 (At |) 2 - pCxj^B^Atj *) 2 

P 2 (t*) = P le Ax-<- + P 2e (l + «5 2 ) + B 5 3 (At*) 2 } - p 3e B 2 (Ax^) 2 

- Pi^jBgAT’r + B 5 (At*) 2 J-F(x 1 )(At*) 2 - F(x^)B 2 (At^) 2 +F(x 2 )At* 

P 3 (t*) = p lg B 5 (AT^) 2 + p^EB^iiT^ + B^Cat*) 2 ] + P 3e { 1_ [ - <£>) 

+ (B 1 + B 2 )](At^) 2 } - p^ e [B 3 AT* + BgB^Ar*) 2 ] 

+ F(x 2 )B 5 (At |) 2 - F(x 3 )At* - F(x^)B 3 (At |^) 2 

P^t*) - P le B i,.(AT*) 2 + P 2e*- B 4 AT * + b 5 (At*) 2 ;1 

+ p^Ax* + P^gCl -[(<» 2 + <*> 2 ) + B_ l ](At^) 2 } -f- F(x 2 )B it (Ax^) 2 

- F(x 3 )(ATg) 2 + fU^At^ (4 8) 


P 5 (-r^) = p 5e [l - B 6 (At^) 2 ] - PggBgAx^ - F(x 5 )at* - F(x 6 )B 6 (ATg) 4 


Pg(T^) = P 5e AT * + Pg e Cl - Bg( Ax*) 2 ] - F(x 5 )(At *) 2 + F(x^)At^ (4 9 ) 

where o>^, a> 2 , At*, At*, B^ and P ie ^ 1 = 1 , ,6, are as defined m 

Appendix E 

In equations (4 8) and (A 9) all of B^, 1 = 1 , ,6, are 

periodic (with a period of 2 jc) and except for B they are nearly constant 

The averages (over the time of one orbit) of B_^, 1 = 1 , , 4 , 6 , were 

substituted into the above equations B was considered constant m 

5 

the time between boundary encounters with its value taken as the value 
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at the time of the previous boundary encounter (its average value was 
not used, since its average over the time of one orbit is zero ) The 
approximate values for p (t*)> i = 1, ,6, which are given by 

equations (4 8) and (4 9) were compared with the "exact" values obtained 
from digital computer solutions 

Equations (4 3) -(4 6) were solved with the aid of the digital 
computer for several sets of values of n , l = 1, , 6 , and several 

initial conditions (some of x(t ) m the interior of S and some 
on the boundary of S) The time of solution was = 6 28 2 ji 

and the solution was printed out at intervals of of 0 01 The motions 

considered were the roll -yaw motions of satellites (l) and (2) and the 
pitch motions of satellite (3) The control strength used was N = 

—It 

2 1 x 10 , = 7 1 x 10 and = 2 06 x 10 These "exact" solu- 
tions for p (t*), l = 1, ,6, were compared with the approximate values 

given by equations (4 8) and (4 $) for as many as ten intervals of time 
per solution as follows Each solution of x (t*), i = 1, ,6, was 

applied through F(x ) 4 o equations (4 8) and (4 9) The initial 
conditions on p^, 1 = 1, ,6, were the same as p_ l (t*)^ i = 1, ,6 

After the initial time t* the values of p , i = 1, ,6, were taken 

to be the values of p (r*), i = 1, . , 6 , at the time the proper (yaw., 
roll or pitch) trajectory projection reentered the proper projection of 
S. These new values of p were retained until the trajectory exited 
and reentered S again 

In most cases the comparison of the approximate solutions with 
the "exact" solutions showed that the approximate solutions differed from 
the exact solutions only m the third significant figure In other cases 

the difference was only in the fourth significant figure Therefore, if 

(4 8) and (4 9) an® substituted into (4 2), the result is a time-varying 
feedback control law which is considered optimal for the present problem 
This control law can he implemented m a controller, but, the devices 

needed to determine if x is exiting S or entering S when |x | = 

-4 1 1 

10 may be complicated and can be unreliable Timers, which are 

activated when x , l = 1, ,6, exit S and reenter S are also 

needed to generate At*, l = 1, ,6, and At* Thus, even if this 
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control law is the optimal minimum fuel control law with hest values of 
n^ and for the above satellites, it might not be the least costly 

control law to implement m practice More is said about this m 
Section B 

Estimates of the "best" values of n and N for the satellites 

1 i 

can be obtained with the aid of (4 8) and (4 9 ) If the control is off 

and x(x*) is exiting S, the time it tahes for |x (x*)[, 1 = 1, ,6, 

-4 -4 1 

to increase from 10 to 1 1 X 10 is generally small enough, so that, 
if HRS^, i=l, ,6, is the no-control right hand side of the equations, 
(43) and (4 5)^ for 'x , 1 = 1, ,6, then 


Ax x pa (MRS )AxJ 


1 — 1, • ,6 


(4 10 ) 


is a satisfactory approximation 

During those periods of time when the control is off, |p (t*)|, 

1 = 2,4,6, must be less then unity For the control to turn on so that 
|x (t*)|, 1=1, ,6, do not increase to values greater than 1 1 X 10 \ 

the values of the corresponding |p (x*)], 1 = 2,4,6, must be unity or 
greater at the time when Ix^x*)] = 1 1 x 10” Since m (4 8) and (4 9) 
p , 1=1, ,6, are nearly periodic with an average over one orbit of 

nearly zero (see Figure 4 lb, for example) a first approximation is 

P 2 (t*) & F(x 2 )At* - F(x 1 )(At ^) 2 


Pj^(t*) «F(x lt )AT^ - F(x 3 )(At *) 2 


Pg(x^) pa F(xg)Axg - F(x 5 )(At |) 2 


(4 11) 


so that for 

which I Ax | 
-4 1 

1 1 X 10 , 


|p (x*) | to grow to unity m the time interval Ax* m 
grows to 10 , 1 e , |x^(x*) | increases from 10 to 

it must be true (as seen from (4 10) and (4 ll)) that 


10 1 = F(x ) w (Ax*)" 1 pa MRS^/lO" 5 , 1 = 2,4,6 

10 1 * F^) pa (AX*)" 2 pa (lffiS^/10" 5 ) 2 , 1 = 1,3,5 (4 12) 
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Equations (4 12) can be used to determine an estimate of the 
"best 1 ' values for n , 1 = 1, . ,6 In (4 12) the largest possible 

values for HRS should be used since smaller values result m n 

which are too small to keep |x | ^ 1.1 X 10 For satellite (2), for 

1 11-4 

example, the largest possible values of HRS for |x | £ 1 1 X 10 

_4 1 1 -4 

are HRS, = 1 1 x 10 , HRS 0 = 1 13 x 10 , HRS~ = 1.1 X 10 , HRS, = 

3 69 x 10 , HRS^ = 1 1 X 10 and HRSg = 2.06 x 10 With these 

values it is found from (4 12 ) that n^ = 2 . 1 , n^ = 1 0 , = 2 1 , 

= 1 6 , = 2.1 and = 3.3 

From several digital computer solutions obtained with the above 

estimates of the best values of n and with several sets of values of 

1 

H , it was found that some of |x |, 1 = 1, . , 6 , grew to values as 

—4 

large as 6 0 x 10 , which, of course, is too large Also, from the 

solutions it was found that for and H^ neither the largest nor 

the smallest of the values tried were best and for H_ the smallest 

0 

value was best. A detailed examination of this result with the aid of 

(4 2), (4 3), (4 5 )) ( 4 . 8 ) and (4 9 ) proved to explain the effects of 

n and H on the solutions and offered new values of n and N to 
11 11 

be tested (This detailed examination is not given here since it is long 
and tedious with many numbers and since similar examinations are given m 
Parts 2 and 3 ) With the newly offered values of and H digital 

computer solutions were obtained It was found that the requirement for 
the "best" values was satisfied by mcreasiing n_^, 1 = 2,4,6, by about 
one over the above estimated values As expected, the accuracy of earth- 
pointing improved and the fuel cost increased to the maximum value, which 
was for no control-off intervals, with increasing n Also, it was 
found that both the fuel cost and the earth-pointing accuracy generally 

increased with increases in H above their "best" values for values of 

1 

H less than about 0 1 For values of N of 0 1 or more the fuel 
1 1 

cost still increased with increases m H but the accuracy of earth- 
pointing dropped off sharply Values of H^ which were smaller than 
the "best" values resulted m the control being on so much longer that 
the fuel cost was generally higher The "best" values of n^ and H^ 
resulted m the pitch axis control staying on about 95$ of the time 
while the roll-yaw control was on about 33$ of the time 
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In summary the main result of this part is the time-varying 
feedback control law given by (4 2) together with (4 8) and (4 9) The 
control law is optimal m the sense that it was obtained from the suffi- 
cient conditions of the maximum principle for the cost functional 



t 1 
o 


and the integral constraint on the states 
rf 

/ f-Jx^T^dT £ A 

T 

O 

Other results of this part are (a) a not-too-expensive method has been 
devised for determining the "best" values of the control law parameters 
and U , (b) the effects of varying the parameters and 

from their "best" values was given, and, (c) satellite (2) has the small- 
est unit fuel cost with the above control Figure (42) shows yaw, roll 
and pitch projections of a state space trajectory of satellite (2) for 
the time of one orbit with the above control and the "best" values of 
n and N 

l l 

The optimal control of this part is not necessarily minimum fuel 
optimal Indeed, it is not difficult to exhibit several controls which 
keep x m S and use less fuel Since it is very difficult and 
expensive computer-time to apply the NCES to the full three -axes satellite 
control problem and since a basis for fuel cost comparison is needed, 
simple single-axis control problems are now investigated 


2 The Steady-State Motion Obtained From PMP for x" = v 


In x" = v let x = 


x 5 , X 


t _ 


x^ and v = Vg Then 


X 3 X 4 


x k ~ V 2 


(4 13) 
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Figure 4 2 Solution of Equations 3)-(4 6 ) for Satellite (2) 
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V 


is equivalent to x" = v For some intervals of time, equations (4 13) 
approximately describe the roll motion of satellites if they are nearly 
symmetric (such that « 0), if they are not significantly affected 
either by the earths atmosphere or the disturbances mentioned m Chapter 
II and if they are controlled such that x = 0, i = 1,2, 5, 6 Of course, 
no such controlling of a satellite 1 s motion is practicable at this time 
However, if v^ is a "coast function" with magnitude N 2 ^ 0 01, equation 
(4 13) adequately describes both the roll and yaw motions during those 
time intervals when the control is In any case, for the purpose of 

comparing the two mathematical methods (approximate and exacts, it is 
worthwhile to study the controlled motion of the simple plant described 
ty (4 13) 

By studying the possible controls and motions which satisfy the 
PCES (see Part 4, Section D, Chapter III), controls which perform satis- 
factorily are obtained and given below (These controls are compared 
later to the controls obtained from PMP with 


J =J [f^x) + MldT 


for the motion described by (4 13), i e , the approximate method ) 

While the phase -plane trajectory is m the interior of S^, the 
roll phase-plane projection of S, the conditions of PMP must be satis- 
fied They are (4 13) and 

P3 = ° 


p 4 = " P 3 
v 2 = h 2 cst(p 4 ) 


(4 14) 
(4 15) 


The parts of the trajectories which coincide with the parts of 
the boundary of S R given by x^ = ± 10”^ are "regular" For a given 
value of N 2 (< «), trajectories which begin m the shaded regions of 

m Figure (4 3) must exit m acquiring equilibrium points, namely 

those points such that x^ = 0 These trajectories do not satisfy the 
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NCES unless the vertical parts of the boundary of S R are replaced by 
segments of parabolic curves (control-on trajectories) The control-on 
parts of the trajectories "A"* TI B" and "C" m Figure (4 3) are such 
curves The trajectories are "regular" All trajectories m Figure (4.3) 
except "B" acquire the part of the line x^ = 0 m or sufficiently near 

In this case the transversality condition (3 l4) with b(x ) 

perpendicular to the x^axis is applied Trajectory "B" is not precisely 

a minimum fuel trajectory except m the limit as e -* 0 However , e ^ 0 

is more practical since imperfections such as tame delays, thresholds, 

etc , (see Flugge-Lotz [12]) are always present m the controller, and, 

for comparison of the NCES with the approximate method, this is an 

interesting case when x^ is free m the interior of S The value of 

N used m Figure 4 3 is the smallest value which results m |x | £ 

1 1 X 10 , i = 3, 4, (a 10$ error m the maximum error) regardless of 

the initial point m £L, The "jump condition" (3 15) is applied at "P" 

Ja 

Since the time required to reach x^(t^) = 0 m S R (or a 
point suitably near so that no future control effort is required for the 
remaining part of the satellite’s lifetime) is not important, the usual 
"bang-cost -bang" control, ie, ± N^ -* 0 -> =F N^, does not result m a 
minimum fuel expenditure (see, for example, Marbach [25] wherein the 
acquisition time is given) 

Now consider the controls and resulting motions obtained from 

PMP with 

T q 

J = J [f-]_(x) + |v[]dT 

T f 

for x" = v and for various times of consideration and compare them with 
the controls and resulting motions obtained from the NCES 

As m Part 1 of this section, the equations of the maximum 
principle for this problem can be written m backwards time as 

' x 3 = " x 4 

= - v (4 1 6) 
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,p 3 ~ 

'p 4 = F(x 4 ) + p 3 (4 17) 

v = H cst(p 3 ) (4 18) 


These equations can be solved m a piecewise manner. 

The piecewise solution of (4 1 6) for intervals of over 

which v(t*) is constant is given by 






(4 19) 


where t * is measured from the last time of control switching, t*, or 

s 

from the initial time, r^, and x 3Q = x 3 (t^), x^q = x^(t*), i = o, s 

If equations (4 17) are solved for the same intervals of r* as 
the piecewise solution of the adjoint variables given m Part 1 and m 
Appendix E, the result is 


P 3 (t*) = P 3e “ F(x 3 )At* 

P 4 (t*) = P 4e + P 3e AT7f ' + f (x 4 )AtJ ~ f (x 3 )(At*) 2 /2 (4 20) 

(Notice that when x e S equations (4 20) are just the backwards time 
solution of (4 14) and that (4 18) and (4 15) are the same for any 
q(x) ) 

With the proper initial conditions equations (4 19) and (4 20) 
can be used together with (4 18) to construct optimal controls and optimal 
trajectories If the boundary of S is not encountered by the trajec- 
tory, which m foreward time goes from a free initial point to a fixed 
final point on the x 3 ~axis ( |x 4 (t^) [ ^ 10 ^), the optimal control 
results m a trajectory which is typically any of those m Figure 4 3 
except the ones which cross the boundary of S This is easily verified 
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with (4 l8)-(4 20), since F^) = 0, p ig = p^t*) = P^I^), i = 3,4, 
and p (t^) is either completely arbitrary (if x(t^) is a given point 
and x(t ) is freely chosen) or such that the transversality condition 
is satisfied (if x(r f ) is any point of x^(-r f ) = 0 and |x 3 (T f ) | £ 10 ) 

If the initial point of the trajectory is so near the parts of 

-4 

by x^ = dfc 10 that an exceedingly large 

control is needed to keep the trajectory m the interior of S_, the 

trajectory is allowed to exit by a suitably small amount As before 

the trajectories are limited by |x | s 1 1 x 10 , i = 3,4 The value 

1 -4 

of used m applying the WOES, namely = 5 0 x 10 , is used 

here 

If the proper final conditions (initial conditions in backwards 
time) are imposed by the adjoint variables and if the "best" value of 

is used, equations (4 l8)~(4 20) give boundary encounter trajectories 
which are precisely the same as trajectory "B" of Figure 4 3 which was 
obtained from the NCES, see Part 4, Section D, Chapter III Figure 4 4 
shows two such phase-plane trajectories Figure 4.5 shows a plot of 
( tt"^) for one of the trajectories as obtained from both the HCES and 
(4 20) The curve m Figure 4 5 which is obtained from (4 20 ) is unique 
only for given p 3 (r^) and p^( T*) Ibis curve is only required to 
pass through the points (t*, +l) and (t*, ~l) Thus, the parameters, 
and ng, can he varied somewhat without affecting the 
control, v^(t*), or the resulting trajectory. The values of these 
parameters used m constructing Figure 4 5 are p 3 (t*) = -0 20, p^(r£) = 
1.004 and n 3 = 0 42 For convenience p^(t*) was chosen to he zero 
so that the "best" value of n^ was the same as obtained from the method 
of estimation of Part 1 

From the above investigation it is apparent that the approximate 
solution, equations (4 l8)-(4 20), gives the same results as the exact 
solution, if it is required that the final point of the trajectory be a 
point of the Xg-axis m S R It is not difficult to see that the two 
solutions also agree precisely m the case that the final point, [x 3 (t f ), 

x 4 

free (m S„) final point requires p 7 (t jj ) = p, (t ) =0 ) 

K or 4 1 


(t^)], is free xn if is exited only once (As m Part 1 the 


p 5 ( ' t S ) ' p 4 (t ? ) 


the boundary of S given 
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Figure 4 4 Nearly Minimum Fuel Boundary Encounter Trajectories for 
x" = V, n 2 = 5 o X 10~ 4 



Figure 4 5 Adjoint Variable m Backwards Time. Note P^(t^) >10 
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If the approximate solution is carried out for t* - t* which 
is sufficiently large for two or more exits of S^, it becomes clear 
that regardless of the values of n^, n^ and N used the approximate 
solution has a fuel cost which is larger than the fuel cost of the exact 
minimum fuel solution The approximate solution can be made better by 

changing the values of and n^ before each exit or S R or by 

applying the solution m a piecewise manner with the jump condition. 

These methods of improving the approximate solutions with more than one 
exit of S R are, generally, not satisfactory since they make the con- 
ditions of the approximate solution as difficult to apply correctly as the 
NCES An alternative to these methods is the use of functions other than 
fnCx) m the integral constraint on the states The functions fg(x) 
and f g (x) given by (3.10 ) and (3 ll) are the same as f-^(:x) when x e S 
so that they result m the same satisfactory solutions if the trajectory 
remains m the interior of S If the trajectory exits the region S 
more than once in the time interval of consideration, the use of fg(x) 
instead of f^(x ) results m a reduction m fuel cost This is true 
since the derivative of f 2 (x) (which appears m the adjoint solution) 
is smaller near the boundary of S than the derivative of f (x) so 
that the fuel cost carries more weight over the trajectory when f^(x) 
is used than when f, (x) is used 

A free final point solution of the approximate problem with 
f^(x) used m the integral constraint is given below The steps leading 
to this solution are very similar to steps leading to the solutions (with 
f^(x) used) given above with the several types of boundary conditions 
For completeness these steps are given m detail 

The equations of the maximum principle for this problem m 
backwards time are 



'x k = - v (k 21) 


I 




r 

10 (x. 


- 


0 



) 


-if. ) 

if x^ > 10 

If |x 3 l £ 10" 4 ^ 


(continued) 
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r 


' P 4 = p 3 + < 




v = H CST(p, ) 


n„ 

10 (*s 

-4 \ 

+ 10 *), 

if 

-4 

x- < - 10 * 

0 J 


n, 

10 4 (x 4 

_4 . 

- 10 4 ), 

if 

x ^ > 10" 4 


0 

; 

if 

U k 1 * 10" 4 

> 

n 4 . 

10 (*4 

+ 10" 4 ); 

if 

1 

S 

V 

(4 22) 





(4 23) 


Equations (4.21) are the same as (4 1 6) and their solution is 
the same as (4 19 ) Equation (4 23) is identical to (4.l8). If the "if" 
parts of the statements which define the functions m brackets m (4 22) 
are assumed but not written, the solution of (4 22) can be written as 


n„ 


r 2 (x 5Q - 10~ if ’)'r if - x^t* 2 + vt* 3 /3 ^ 


P 3 h*> - - 4 


0 


2(x 3Q + 10~^)r* - x^qT* 2 + vt * 3 /3 


} + P 3e 24) 


J 


r 


n„ 


P^*) = “ 


( 0 


4. 


+ P 4e + p 3e T * 


n. 


+ 10 t* < 0 


: 30 " 

10“ 4 )t* 2 

1 

** 

O 

'30 + 

io-^t* 2 

‘ x 4o 

T* 



' x 4o 

- 10 - 

vrt /2 

0 



x 40 

+ 10* 4 - 

vt*/2 


J 


(4 25) 


where t*, x^q and x^ are initialized at each time of encounter with 
the boundary as well as at each time a switch m control occurs 

Suppose that the initial point is given by [x„(r v ), x, (t*)] = 

j- O O 4 0 

(-5 X 10 10 and m forward time is a free final point In this 
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case p^ e = P 3 (t*) = 0 and = p^(t*) = 0 Until the boundary of 

S is encountered (point "1" in Figure 4 6s.), sat at t*, the solutions 
(4 19), (4 23 ), ( 4 . 24 ) and ( 4 . 25 ) give 

x 3 (t*) = x 30 - x^t* - 5 X 10“5 _ io" 5 t* 

V T * } = x l*o = 10 “ 5 

v = 0 

P 3 (t*) = 0 

p ^(t*) = 0 (4.26) 

From the first equation of (4 2.6) the time of encounter with 
x 3 (t*) - - 10" 4 is found to be = 5 Since x 3 ,x^,p 3 and p^ are 
continuous, they have the same values at ( T q) + as (t£) Thus, 
until the control turns on, say at t*, the solution is 

x 3 (t*) = x 3 q - = -10"^ - io“7* 

x^(t*) = x^ Q = 10“ 5 

f 

v(t *) = 0 

P 3 (t*) = - 10 72{-x 4o t^3 = 10 5 t* /2 

p^ 7 *) = - 10 3 /aC-x^-r* 3 /4) s 10 3 t * 3 /8 (4.27) 
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With N = 5 x 10 the allowable error will not he exceeded, 

I I «!}. 

i.e , ^ 1 1 X 10 , if the control turns on when x„ = - 1 08 x 

_i|. 6 3 _i, 

10 From the first equation of (4 27) with x_(t*) = - 1 08 x 10 it 

O ^ 

is found that t* = 0 8 (measured from t*) For the control to turn 
on at r* } must grow to unity at r*. Thus, from the last equation 

of (4.27) the "best" value of n^ is found to be n^ = 6 2 From the 
fourth equation of (4 27) it is found that p (t^) = 5-0 Thus, for the 
next part (from point "2" m Figure 4 6a), of the calculation x 3Q = 

-1 08 x 10“\ x = 10~ 5 , P 3e = 5 0, p 4e = 1 0 and v = TS = 5 X 10 -4 

Until the boundary of S is again encountered (point "3" m Figure 4- 6a) , 
say at t*, the solution is 

x 3 (t^) = Ut * 2 - + x 3Q = 5 x 10“\*- 2 + 10 ~^t* - 1 08 x 10 ^ 

x^(t*) = - Wt* + x^ 0 = - 5 x 10 + 10 ^ 

v(t*) = W 

P 3 (t*) = - 10 3 / 2 {2 (x 30 + 10"V - x^r* 2 + vt* 2 /2} + p 3e 
- 5 0 + 12 5r* + 7 8 t * 2 - 130. Or* 5 

Pj^(t^) - - 10 3 /2 £(x 30 + 10 _i| ')T^ 2 - x^ q T^ 3 /3 + vt*^/ 12] + Plfe + p 3e r* 

= 10+5 Or* + 6 25t + 2 + 2 6 t * 3 - 32 5-r*^ (4 28) 

The time of reentry is found from the first equation of (4 28) 
to be t* = 0 126, since x 3 (t^) = -10" : From (4.28) with t* = 0 126, 

the initial values for the next part of the calculation are found to be 
x^ 0 = - 5 30 x 10~^, p 3e = 6 43, p^ e = 1 73 and v = N The next part 
of the solution is such that p^(t*) continues to increase so that the 
control stays on and drives the trajectory to the boundary again Thus, 
until the boundary is encountered again (point "4" m Figure 4 6a) , say 
at the solution is 
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x 3 (t:*) = Nr* 2 - + x 3Q 

_4 2 -5 -4 

= 5 X 10 x* + 5 3 X 10 - 10 

x^(t*) = - m* + x^ 0 = - 5 X 10“^T* - 5 3 X io~ 5 

v = S 

P 3 ( ^ } = P 3e = 6 

p ^(tt*) = p^ e + P 3 e 'f* = l 73 + 6 43t+ (4 29 ) 

The time of exit, 
x^ = 0 094 since X ^( T ^) = 
initial values for the next part are found to be x^q = -9 1 x 10 ^ 

P 3e = 6 43, p^ e = 2 33 and v = N The next part of the solution con- 
tinues until the control switches, say at t*, and is 

X (t*) = 5 X lO-^T* 2 + 10"^T* - 9 1 X 10“ 5 

O 

X. (t*) - - 5 X 10 T* - 10- 4 

P 3 (T *) = P 3e 

n, . 

Pi,.(t*) = P^g + P 3e T * + 10 t ^ x 4o + 10 ~ “ VT */2} 

n, 1 

= 2 33 + 6.43T* + 10 t*{- 5 x 10 #/2) (4 30) 

The control switches to zero when p^(t*) =10 and the trajectory 

becomes parallel to the x,-axxs Since it is required that |x. (t*) | £ 

-4 5 ^ 

1.1 x 10 , the control switch must occur no later than r* = 0 02 Thus, 

from the last equation of (4.30), it is found that I 1 is required 

If n. = 7 1 and t* = 0 02 are used m (4 30 ), the initial values for 
45 _ 5 -4 

the next step are found to be x 3Q = -8 88 x 10 , x^q = -1 1 X 10 , 

P 3e = 6 43, p^ =10 and v = 0 The solution for the next part (from 

the point "5" m Figure 4 6a), is 


t|£, is found from the second of (,4 29 ) to be 
-10 ^ From (4 29) with xff = 0 094 the 
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x 3 (t*) = x 30 - x ^.o T * = -8 80 X 10 5 + l l x io""\* 
x^Ct*) = x^ 0 = -1 1 x 10"^ 


P 3 (t*) = p 3e = 64 3 

n, . 

p^r^) = p^ e + p r* + 10 t*{x^ q + 10 -4 } =10- 118 6t* (4.31) 

This solution is for the time interval (t* - , x£] where xg is the time 
when the control switches from zero to V = - N This switch m the con- 
trol occurs when p^(rg) = -1 Thus , from the last equation of (4 3l) it 

is found, that = 0 017 The initial values for the next step can he 
found from (4 3l) with xp = 0 017 They are x„~ = -8 69 X 10 - ^, 

_k ° 5U 

x^ 0 = -1.1 x 10 , = 6.43, p^ e = -10 and v = - N Thus, the 

solution for the next step, which is until r* at which time the trajec- 
tory reenters S, is 

x 3 (t*) = - Hr* 2 - x^t* + x 3Q = - 5 X 10 -if T* 2 + 1 1 x lO^T* 

- 8 69 x 10" 5 

x^(t*) = Nt* + x^ 0 = 9 x 10 - \* - 1 1 x 10"^ 

P 3 (t * ) = P 3e = 6A3 


P 4 M = P 4e + P 3 e -r* + 10 4 t*{x 4o + 10‘ 




1 0 - Il8r* + 3120 t*^ 


(4.32) 


The time x* can he found from the second equation of (4 32) 
since x^(t^) = - 10 It is t* = 0 02 The initial values of the 
next (and final) step can he found, from (4 32) with x* = 0 02 They are 

X 30 = - 8 4 x 10 “ 5 , x 4Q = - 10 "\ p 3e = 6 43, p^ e = - 2 12 and v = - U 
The solution for the last step which is until rg is 

-5 


x 3 (t*) = - 5 x 10 ^t* 2 + io“^t* - 8.4 x 10’ 


(continued) 


70 



(4 33) 


x^Ct*) = 5 X 10“*V - 10“ 4 
P 3 (t*) =6 43 
P^(t*) = - 2.12 + 6 43t* 

The final time T f = T * + • + is chosen so that the trajec- 
tory begins (in forward time) m S The solution is given m graphical 
form m Figure 4 6 

In summary the main results of this part are (l) minimum fuel 
controls for the single-axis motion of satellites if the motion is 
described by x = v, (2) the demonstration of the fact that the solutions 
obtained from PMP with 


J = 



(f(x) + j V | )&T 


T 

O 


(f(x) = f-^(x), ^2^—^ ’ } axe m - Lnimum solutions if the 

boundary of S is exited no more than once, (3) the use of smooth 
functions like fg(x) results m less control effort than the piecewise 
smooth function f^(x) if more than one exit of S occurs 

p 

3 The Steady-State Motion Obtained from PMP for x" + a x = v, 
a 2 > 0 

2 p 

The equation x" + a x = v describes a stable system if a > 0 
If the pitch motion is controlled so that Ix^t)] S' 1 1 x 10”\ 1 = 5,6, 
and if the satellite has inertia properties such that k^ 1 and k -1 
(like satellite (2), for example), then both the roll motion and the yaw 

O P 

motion are very nearly described by x + a x = v with a > 0 If the 

-6 

eccentricity of orbit is less than 5 X 10 - , k„ 1 and the yaw-roll 

O 1 

motion is controlled so that |x (t) | sll x 10 , 1 = 1, ,4, then the 

1 p 

pitch motion is also very nearly described by x" + a x = v Thus, the 
application of the theory of Chapter III to 


x‘ = x,. 


xl = - 


2 

a X-, 


+ v 
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leads to useful information about the theory and the required control 
In this part of Section A, the requirements on the control and the 
methods of analysis are the same as in Part 2 of Section A 

First consider the necessary conditions for exact solutions 
(NCES) of Part 4, Section D of Chapter III While the state-space trajec- 
tory is in the interior of S the conditions of PMP, i.e , (4 34) and 



P' = - p x h 55) 

v = 3BT CSl(p 2 ) (4.36) 

must be satisfied -with the proper boundary conditions 

Three representatives minimum fuel phase plane trajectories which 
do not encounter the boundary are shown m Figure 4 Ja. These trajectories, 
denoted by A, B and C, are the results of applying PMP with three different 
sets of initial conditions on (xy,Xg) and (Pq-’Pg^ ' Hie circular region, 
r, which is just inside the yaw projection of S, Sy, is the region to 
be acquired if initially (after acquisition to S) the trajectory is not 
already m T The region r was chosen to be acquired since once it is 
acquired no future control effort is needed to keep the trajectory m Sy 
(in the absence of detrimental disturbances) Also, less fuel is needed 
to acquire P (from within Sy) than any other circular region m Sy 
The final points, P, of the trajectories, A, B and C, correspond to the 
final time of consideration, namely t , For trajectory A the final point 
is fixed The final point of trajectory B is free so that Pq( T f) - 
p (t ) = 0. The trajectory C acquires the boundary of r so that the 
transversality condition applies 

There are points m Sy from which P cannot be acquired without 
encountering the boundary of Sy A typical case is shown m Figure 4 8 
The final point of the trajectory is fixed so that p^ ( ) , i = 1,2, are 
entirely unknown a priori. The jump condition is applied at the junction 
point denoted by the number "3" (Those parts of the boundary of Sy 
which are vertical lines are replaced by arcs of circles about (:Htf,o).) 

Xn Figure 4 8a the normal to the boundary at the junction point makes an 
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angle of 3° with the x-^-axiS; so that; m Figure b 8b the "jump" m 
(P^P 2 ) at "3" is at an angle of 3 with the p^-axis In this appli- 
cation the jump condition constant is negative. The normals to the 
vertical boundary lines at the boundary points "1" and "2" are parallel 
to the x^-axis. 

The part of the trajectory m Figure 4 8a which begins at the 
junction point "3" and ends at the fixed final point P is part of the 
optimal trajectory from "3" to the origin for acquisition times greater 
than 5it/8 (See Marbach [25]; Fig 2 11.) This part of the trajectory 
is minimum-fuel optimal since it is part of an optimal trajectory The 
time of acquisition (> 5tc/8) is taken as large as possible since the 
fuel cost decreases with increase m acquisition time (See Marbach’ s 
Fig 2 15 ) 

2 

For x" + a x = v as for x" = v m Part 2 it can easily be 
seen that the optimal (minimum fuel) controls obtained from the NCES 
with the jump conditions applied at boundary encounter points are the 
same as the optimal controls obtained from PMP with 


J = 



[f(x) + |v|]dr 


T 


o 


The equations of the maximum principle for 


J = 



[f (x) + I'vjJdT 


T 

O 


m backwards time are 


i 


x. 



2 

a X 1 ” T 


'P-L = - F(x x ) - a 2 p 2 
’p 2 = F(x 2 ) + p 2 


(4 37 ) 


(4 38) 
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(4.39) 


v = 3Sf CST(p 2 ) 

As m forward time the solution of (4 37) is a parametric representation 
of the equation of a circle m the (x^x^) -plane The center of the 
circle on the x^-axis at = 7 (t = ± I or zero) The trajectory- 
proceeds m a counter-clockwise sense with increasing t* The piecewise 
solution of (4 38) can he obtained by the method of Appendix E and is 

Pi(t^) = p cos(aAir*) - P 2e a sin(aAT*) - P(x 1 )a sin(aAt|) 

- F(x 2 )[1 - cos(aAT*)] 

= ? 2 e cos ( a 4' 1 "^) + Pp e a ^ sm(aAT*) + F(xg)a ^ sin(aAT|) 

- F(x 1 )a -2 [1 - cos(aAT^)] (4 40) 

Figure 4 9 shows a graph of a solution of (4 38) as obtained 
from (4 40) This solution corresponds to the solution of (4 37) shown 
m Figure 4 8a m forward time (as indicated by the arrows) As before 
m Part 1 and Part 2 of this section., the values of n^ and n 2 depend 
on the initial values of the adjoint variables. The initial values of 
p^ and p 2 used m the solution shown m Figure 4 9 are a worst case 
choice Although Figure 4 9 shows a worst case, there is still some 
resemblance to the curve of Figure 4 8b, which is obtained from the KCE8, 
if differences m scaling are accounted for If the initial (backward 
time) value of p^ is made more negative, the solution (4 40 ) of the 
adjoint equations which corresponds to the trajectory of (4 8a) agrees 
perfectly with the curve m Figure 4 8b for all except about one -fourth 
of the time of consideration The difference m the curves is due to 
the mechanism which causes the change m the adjoint variables on encount- 
ering a boundary For conparison purposes Figure 4 10 shows a solution 
obtained from (4 40) which mostly agrees with the "jump" solution 

The above example, which show excellent agreement m the solutions 
obtained from the NCES and from PMP with 
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are representative of all optimal solutions which must (for a given value 
of h) exit the region Sy (by more than some small allowable amount) once 
and only once The mam reason this generalization can be made is the 
fact that the parameters n^ , n^ and ji are freely chosen 

O 

4. The Steady-State Motion Obtained From PMP for x" + a x = v, 
a 2 < 0 

2 2 

The equation x" + a x = v with a <0 gives an accurate 

description of the yaw motion of satellite (3) for practical time intervals 

if the roll and pitch motions are suitably controlled This equation also 

describes the pitch motion of satellite (2) near apogee and perigee if 

the yaw and roll motions are suitably controlled 

The equations of the maximum principle for this part are the 

same as m Part 3 except that here a <0 [See (4 34)-(4 36) ] Figure 

4.11 shows some phase plane plots of solutions of (4 34) and (4 35) for 
2 

a = -1 The curves are hyperbolic 

Ideally, i.e , m the absence of disturbances and imperfections 
m the controller, the best control is the minimum- fuel control which 
acquires the line PQ m Figure 4 11a (if the line PQ is acquired m 
this ideal case, no future control effort is required, since once PQ is 
acquired the trajectory moves very slowly along PQ to zero ) Since the 
transversality condition must hold at the final time, which is free, the 
minimum- fuel control which acquires the line PQ must be on ( Htf ) all of 
the time of acquisition (See Figure 4 lib ) Therefore, m this ideal 
case the control must be -N above PQ and must be +F below PQ 

The magnitude of the control used m constructing Figure 4 11a 
is too small For two reasons N should be larger They are (l) the 

error m x n ( t) grows to undesirably large values, i e , lx (t) | > 

_4 1 -4 1 -4 

1.1 X 10 , if initially both x^. and x g are near 10 or -10 , 

(2) for most of the cost of acquiring the line PQ decreases with 

increase m N until the cost is the same as for x" = v m acquiring 

the same line This last reason is easily validated by considering the 






time integral of the last of the system equations (equations (4 34)) 

± = x 1 At ± NAt (4 4l) 

where is the mean value (of the integral mean value theorem) for the 

time interval At, Ax p > 0 and N > 0 

X 2 


Xi 


Figure 4 12 Regions of Increase m Cost With Increase m N 

If (4 4l) is solved for IN At, whose magnitude is the cost, the res ult 
is 

U = INAt = ±Ax 2 - x^At (4 42 ) 

As N increases the value of x. approaches the initial value of x, 

/ 1 -4 1 

(Even for small N, say N = 2 x 10 , the change m x^ is small 

compared to Ax g for most initial conditions and time intervals See 

Figure 4 11a ) Since the time interval required for x^ to change by a 

given amount decreases with increase m N, it is easily seen from (4 42) 

that the cost decreases with increase m N whenever the change m Xp 

is positive or negative) and x ± is negative (or positive) In the 
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limit J = AXg In the shaded regions of Figure 4 12 the cost is less 
than AXp for bounded Ef and m the limit J = Ax^ 

In the more realistic case m which imperfections m the con- 
troller and disturbances of the motion are present, the line PQ, cannot 
be precisely obtained, and, even if it could be acquired, the control 
effort is not exactly a minimum The disturbances of the motion are 
mostly due to the coupling with the motion about other axes, the periodic 
oblateness part of the earth’s gravitational torque and (rarely appearing) 

meteoroids These disturbances can be small enough for the trajectories, 

-6 

which are never nearer than about 10 to the diagonal no-disturbances 

trajectories, to remain nearly hyperbolic In this case the control law 

described by Figure 4 13a is very conservative with respect to cost 

(With 4 x 10 ^ the control is on only about one-seventh of the time 

of one orbit ) The trajectory m Figure 4 13a is realistic m the sense 

that a time delay of 0 1 sec is assumed m the controller and a small 
-6 

amplitude («s 10 ) sinusoidal disturbance is assumed to affect the motion 

of the system described by equations (4 34) with a = - 1 The cost of 
the limit cycle motion of Figure 4 13b is about ten times smaller than 
the cost of the controlled motion of Figure 4 13a However, the limit 
cycle motion requires sensors (of the state variables) which give one 
order of magnitude greater accuracy If the two triangular control-on 
regions m Sy which contain part of the x^-axis are replaced with 
control-off regions, the control effort is not significantly increased 

unless a large, but rare, disturbance causes the trajectory to acquire a 

-4 

point near x^ = 0, x^ = ± 10 This simplification of the control 
logic can result m a lower overall controller cost even 'though it can 
cause an increase m control effort (See Section B) 

As m Part 2 and Part 3, the solution of the optimal problem with 

r f 

J = J [fy(x) + |vf]dr 

T 

o 

is precisely the same as the solution obtained from the NCES except for 
the discrepancies m the adjoint variables which occur during the time 
when the state trajectory is outside of The adjoint variables for 
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a. A Low Cost Control Law and Trajectory for a Realistic Case 



b A Low Cost Control Law which Results m Limit Cycle Motion 

Figure 4 13 Low Cost Control Laws and Trajectories., N = 0 05, Time 
Delay =01 sec 



this problem with 


J = J [f-Jx) + |vf]dT 

T 

O 

can be obtained as functions of t* as m Appendix E and are given by 

*i 

P-j^T*) = p le cosh(cAr*) + P 2e c sinh(cAT*) - F(x^)c“ sinh(cAT*) 

+ F(x 2 )[cosh(cAT*)-l] 

1 "I 

P 2 (^) = p 2e cosh(cAT^) + P le c sinhCcAr*) + F(x 2 )c sinh(cAT*) 

- F(x^)c 2 [cosh(cAT*)-l] (4 43) 

where c = - a >0 (it should be noticed that equations (4 43) are 
the "unstable" counterparts of equations (4 9) which are for the "stable" 
pitch motion ) 


B. THE STATION-KEEPING CONTROLS - GENERAL ORBITS 

1 Station-Keeping Control of Satellites m a Maximum Gravitational 

Torque Orbit 

In Section A the orbits were restricted to nearly polar or nearly 
equatorial orbits which were of sufficient altitude that the aerodynamic 
torque was insignificant. 

If the orbits are not nearly polar or nearly equatorial, the 
oblateness terms m the gravitational torque can be very significant m 
the case that high-accuracy earth-pointing is needed. Consider equations 
(2.17) with the parameters of satellites (2) and (3), for example. For 
all but about one -tenth of an orbit the terms which contain either c. 


'0 


-4 


or Sq as factors are larger than the other term when jx | ^ 1 1 X 10 , 

i=l, ,6 Thus, as a limiting case the equations of steady-state 
motion for satellites such as (2) and (3) are 
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x 2 = “ 2k l A 3 J(r E /r)2s 6 C S C e + V 1 
X 3 = X 4 


x 4 = 8k 2 A 3 j(r E /r)2s & C 8 S 0 + V 2 

x 5 = x 6 

= -A 2 + v 3 (if. 44) 


Since e - 0 
so that for 



01 m these cases, and (r^/r) are very nearly constant 

8 = 45° equations (4 44) are closely approximated by 


cos(t +0 ) + v^ 


sm(r h 0 ) + Vg 


sm(r + 0 q ) + (4 45) 


•z rz 

■where for satellite (2) = -1 45 x 10 , Eg = -5 05 x 10" and for 

satellite (3) E ± = 1 26 X 10~ 5 , Eg = 5 22 X 10 -3 

In contradistinction to the approximate motions of Section A, 
no equilibrium point exists an the motion described by any one of these 
three pairs of equations (which m this limiting case describe the yaw, 
roll and pitch motions) In this case there exists no region of stability 
m S (which can be acquired with the aid of a minimum fuel control) 
which is such that no control effort is required to keep the trajectory 
m S for the remaining lifetime of the satellite 

Except for the equations of motion., which have a sine forcing 
term, the equations of the maximum principle for each pair (4 45) are 
the same as for x" = v The adjoint equations m the case of an integral 
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constraint on the state variables are (m backwards time) the same as 

(4 18) 

Although the sine forcing terms are no longer than 01, their 

i i “4 

effects are most important when |x | s 10 so that the optimal feed- 
back station-keeping control law must be highly time varying This fact 
and the fact that no equilibrium points and no region of stability m S 
exist make the determination of practical, i e , suboptimal minimum-fuel 
control law improbable However, enough of the characteristics of the 
optimal control law are known so that a control which performs satisfac- 
torily can be devised 

The optimal control must be a "coast function" as given by (4 19 ) 
If the time and the states at which the control optimally switches can be 
found, the optimal feedback control will be known If the NCES are applied 
to the problem, two characteristics of the minimum fuel control become 
apparent They are (l) the initial (and final) values of the adjoint 
variables, as well as the jump in the adjoint variables at junction 
points, must be small (about one or less m and ten or less m ) 
m most applications for low fuel cost and for minimum dispersion of the 
trajectory, (2) the switches m low cost controls occur when the trajec- 
tory is near the boundary of the station-keeping region 

In the maximum principle is applied to the minimum-fuel problem 
of station-keeping via an integral constraint on the state or equivalently 
to the optimal problem with 

T f 

J = j [fqOk) + |v|]dT 
T 

O 

it is found that for the "best" values of n and N the switches m 

l i 

the control occur when the trajectory is near the boundary However, for 

i i -4 

the values of n^ which keep |x | £ 1 1 x 10 , the pitch control, 

Vg, is on about 95$ of the time of consideration (The time of consid- 
eration was 2jc, the time of one orbit ) In this limiting case the 
controls of the yaw and roll motions are on about 90$ of the time when 
the "best" values of n^ and IT are used (By reducing the values of 
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n and allowing more error m earth-point ing the fuel cost can he 
reduced ) 

Since the solutions obtained from the optimal problem with an 
integral constraint and from the NCES of the minimum fuel problem all- 
exhibit characteristic (2) given above, it is logical to ask: should 

successive switches m the control occur when the trajectory is near one 
edge of the boundary or should successive switches occur when the trajec- 
tory is alternately near two edges of the boundary 9 (in the case of the 
NCES either is possible, but, m the case of the integral constraint the 
latter is generally the case ) The analysis which follows shows that 
both methods of control switching result m the same fuel cost over the 
lifetime of the satellite 

Take, for example, the third pair of equations (4 45) The 
solution for xg(r) is 

Xg(r) = Xg Q + 2e(cos @ Q - cos(t + 0 Q )) + v 3 (t - t q ) (4 46) 

where t ,0 ,X/-„ are initialized after each switch in the control Let 
o o 60 

At = t - t q and Axg = x^(t) - xg Q . Then, since At ^ 0 15 (and is 
usually much less), an excellent approximation of (4 46) is 

±Axg (v 3 2e Bin 0 q )At (4 47) 


where the upper signs correspond to 0 ^ 0 q < it and the lower signs 
correspond to it s < 2 ji If the control switches consistantly between 
off and on when xAr) is one of two given values (eg, both of which 

0 -4 -4 

are near x r = 1 0 x 10 or x,- = -1 0 x 10 or one of which is near 

6 -4 6 _4. 

Xg = 1 0 x 10 while the other is near Xg = -1 0 x 10 ), then the 
ratio of the time the control is on to the total time of consideration is 
found from (4 47) to be 


At 


ON 


At on + At q FF 


It - sin 0 

N 3 o 


(4 48) 


Thus, for a given control magnitude, N 3 , the fuel cost depends only on 
9 q and e and is independent of the values of xg at which control 
switching occurs so that successive switches which occur when the 
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trajectory is near one edge is no more preferable than successive switches 
at alternate edges 

A similar analysis of the dependence of the optimal switching of 

the control on x K was made For any given value of xs, \x,- | £ 10 \ 

| | -4 o o 

no values of x<- } |x^ | ^ 10 , at which optimal switching should occur 

were found The reasons for this are considered to be* (l) the rapid 
changes m the state due to the forcing term., (2) the absence of a pre- 
ferred point or region to acquire , (3) the periodicity of the forcing 

terms 

In summary it has been found for the limiting case of this part 
of Section B that at least one of any two consecutive switches m the 
minimum-fuel control must occur on the boundary of S and there are no 
preferred sets of values of (x,x‘ ) = (x-^x^) or (x 3 ,x^) or (x,-,xg) 
at which the other switch m the minimum-fuel control should occur Thus, 
for the motion described by any one pair of the three pair of equations 
(4.45) the control law given m Figure 4 14 is expected to result m 
satisfactory steady-state performance The distance d / 0 is not a 
consequence of the theory but is an innovation which is made necessary 
by the very imperfection m the controller which makes the control law 
workable — the control time delay For a given d^O (d £ 10 ) and 

a certain size time delay (or control magnitude) there is a maximum value 
of control magnitude (or time delay) at which chatter with no control-off 
intervals will occur If the control magnitude (or time delay) is greater 
than this maximum value for the given value of d j 0, fuel is wasted 
(See Part 4 of this section ) 

This control law, although simple, also results m satisfactory 
performance (see Part 5 and Chapter VI ) when the periodic exponential 
forcing function of the aerodynamic torque is dominant as m the ultra 
high-accuracy earth-pointing of satellite (4) (in this ease, as above 
m the case of the sine forcing term, it is not difficult to show that 
there are no preferred states m S at which control switches should 
occur ) 
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2 Improvement of the Control Law of Part 1, Section A 

Although the suboptimal station-keeping control law., equations 

(4.2) with the second and fourth of equations (4 8) and the second of 

equations (4 9 ), was derived for "stable " satellites m either polar or 

equitonal orbits., similar control laws can be derived for most "unstable" 

satellites m general orbits All of these control laws are similar m 

the terms m which the derivatives of the "penalty function" , F(x^), 

i = 1 ; , 6, enter These terms are zero when x is m S All 

other terms contain the piecewise constant variables; p ; 1 = 1, 6, 

10 

whose average values over several orbits are nearly zero 

If p ,i=l; ; 6 ; are taken as identically zero; the 

10 

control law is greatly simplified, and, for the "best" values of n 
and IT it is, m each phase plane projection, very nearly the same as 
m Figure 4 lb with d = 0 In comparison, the fuel cost m the case 
p = 0 was found to be 25% (roll -yaw) to 66% (pitch) less than when 
both the full system equations and the full control law was used The 
largest percentage values correspond to those satellite configurations 
which require the greatest control effort for earth pointing 

3 The Controls — Approximate Motions 

The control law for single-axis motion which is given m Figure 

4.1b is expected to result m the efficient control of the approximate 

motions of this section and Section A of this chapter Since the most 

general nearly earth-pointing motions of the satellites considered here 

can be approximated (in a least a piecewise sense) by the approximate 

motions above in. this section and m Section A, the three-axes control 

law whose components are as m Figure b lb is expected to result m the 

efficient station-keeping control of the general motion 

If the control law given m Figure b lb is used m the control 
2 2 

of the motion of x" + a x = v, a >0, and if a realistic time delay 
m the control is assumed, say 0 1 sec , the trajectories m Figure b 15 
are typical These trajectories compare well with the optimal trajectories 

of Part 3, Section A (See Figure b.7) ) The distance, d, was taken 

-b 

to be zero, the worst value If d = 1.0 x 10 , the agreement with the 
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optimal trajectories of Part 3, Section A is even better and the fuel 

cost is reduced m some cases by as much as kj$> However, if no future 

control effort is needed, the fuel expenditure is minute even with d = 0 

2 2 

Consider the motion of x + a x = v, a <0, when v is as 

represented m Figure 4.l4 If the time delay in the control is 0 1 sec 

2 

the trajectories m Figure k l6 are typical m the case a = - O.95 

(This value of a is slightly more realistic for approximating the 

"unstable" yaw and pitch motions than a = -1 ) The distance, d, is 

taken just large enough so that the control (with a time delay of 0 1 sec 

is not on all of the time if chatter motion near the two corners (x,x' ) 

, -k -lu 

(HO , HO ) occurs For minimum fuel cost d should be as small as 

-5 

possible The fuel cost m this case with d = 10 and N = 0 05 is 
about 

r f 1 1 -4 

J = / IvjdT = 3 X 10 

T 

O 


per orbit 

b The Controls — Station-Keeping Motions of the Satellites 

The four satellites as described by equations (2.17) or 
equations (2 18) with the control for each axis as described m Figure 
4 l4 were simulated on a PACE TR-48 analog computer (The patching 
diagrams are given m Appendix F ) Some of the results of the many 
simulation runs are shown m Figure b 17 — Figure k 20 as phase plane 
trajectories (or more precisely, the projection of trajectories into the 
phase planes) (The plots were made by an x-y plotter which as directly 
connected to the computer output ) In all figures the initial point, 

■f 

P , is outside of S and on the boundary of a region, S , which is 

twice the size of S (m maximum dimension). The reason P is not 

a + 

taken m S but on the boundary of a newly defined region S is given 
m the acquisition chapter. Chapter V The real time (satellite time) 
of each run is at least the time of one orbit When more information is 
obtained from a trajectory corresponding to more than one orbit (for 
example, the stabilization of the roll-yaw motion of satellite (2)), the 
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Figure 4 1 6 The Control of Figure 4 14 Applied to x" + a 2 x = v, 

a 2 = - 0 95, u = 0.05, d = 10~ 5 
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2xl0' 4 X 5 


T- - t = if-Jt rad 
f 1 

(two orbits) 

9 = jt/2 rad 

o ' 


c Pitch Trajectory, = 0 2, t~ = 0 05 sec , t^ = 0 01 sec., 
d = 2 5 X 10“5 


Figure 4 17 Continued. 


9k 




+ ^ 




X 6 



t_ - t = 2ic rad 
f 1 

(one orbit) 

0 = 90° 

o 

c. Pitch Prajectory., U_ = 0.03., t~ = 0 2 sec , t* = 0 04 sec , 
c o a a 

d = 2 5 X 10"- 7 

Figure 4 l8 Continued. 
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t_ ** t = 2jc rad. 
f i 

(one orbit) 

0 q = jc/ 3 rad 

c Pitch Itajectory., = 0.03, t^ = 0 2 sec , t* = 0 Oit sec , 
d = 2 5 X 10“ 5 

Figure 4 19 Continued 
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T_ - t = 2it rad 
f 1 

(one orbit) 

9 = 0 rad 

o 

c Pitch Trajectory, U 3 = 0 12, t"=02 sec , t* = 0 Ok sec , 

a = 2 5 x 10 “ 5 

Figure 4 20 Continued 
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real time is the time of two orbits The point m the orbit correspond- 
ing to the initial time is given by 0 qJ the angle of the satellite from 
perigee The time delay m each part of the controller is given m each 
figure as t^, which corresponds to -U, and t*, which corresponds to 
+N It should be noted m reading the figures that the lightest regions 
of the rapid chatter parts of the pitch trajectories generally correspond 
to the most rapid chatter while the darker parts of the trajectory corre- 
spond to slower motion. In Figure 4 20a is shown a yaw trajectory of 
satellite (4) (which is greatly influenced by the aerodynamic torque) for 
the time of two orbits. The part of the trajectory corresponding to the 
second orbit overlaps part of the first orbit trajectory However , the 

early rapid chatter motion, which occurs when the aerodynamic torque is 
near its maximum magnitude, is light enough for the later less rapid chatter 
motion to be seen over it. In the roll and yaw parts of Figures 4 17 — 

4 20 the boundary lines of the phase plane projections of S are not per- 
fectly straight as m Figure 4 14 These boundaries were created m the 
analog simulation by an electronic signum function generator which was not 
perfect, 1 e , the characteristics of the diodes used were only nearly 
ideal The large control magnitudes required m pitch resulted m chatter 
motion along the nearly vertical parts of the boundary unless these parts 
of the boundary were straightened up somewhat Therefore, for the pitch 
parts of Figure 4 17 — Figure 4 20 the boundary lines were (electronically) 
made very straight by connecting two diode signum function generators m 
series 

Table 4 1 gives the fuel expenditure (nondimensional) for the 
simulation runs of Figure 4 17 — Figure 4 20 and Jp are the 

total fuel expenditures for the roll-yaw and pitch controls, respectively 
tL, A _ A , for example, is the cost of acquiring the roll-yaw projection 
of S from the point P JL oa , for example, is the cost of the 
steady-state pitch motion Hie fuel expenditure for these runs were 
typical of all runs of the same duration The average (the initial points 
m S and the angle were varied) steady-state fuel expenditure for 

one orbit differed from the steady-state values m Table 4 1 by about 
two or three percent Of course, the acquisition part of the fuel 
expenditure depended greatly on P and to a lesser extent on 
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TABLE 4 1 FUEL EXPENDITURE FOR THE SIMULATION RUNS CF FIGURE 4 1? - 
FIGURE 4 20 


SATELLITE 

(1) 

(2) 

(3) 

00 

J R-Y 

# 

7 52xl0" 3 

* 

7 60 xlo “ 3 

3.71X10 " 3 

* 

7 88xl0 -3 

J P 

* 

3.05X10 " 1 

6 03xl0“ 2 

7 52X10" 2 

* 

13 94xio" 2 

J R-Y*ACQ 

2 88x10" 3 

4 34xl0“ 3 

1 25 x 10** 3 

2 04xl0" 3 

J P*ACQ 

3 09X10 " 2 

-4 

3 0x10 

1 51X10 " 2 

1 03X10" 2 

(t) 

J R-Y* SS 

-f 

4 64x10" 3 

■Sf 

3 26xl0 _ 3 

-3 

2.46x10 

5 84x10" 3 

(t) 

J p* SS 

* 

2. 74x10 “ 1 

6 03X10" 2 

6 oixio" 2 

* 

12 91X10 " 2 


■^Denotes two orbits* otherwise one orbit 

NOTE The nondimensional cost values are translated m Chapter VI into 
the number of pounds per orbit (or year) for the example satellites of 
Chapter II 

( t ) RECALL • Satellite (l) is "mildly stable" m all axes but the pitch 
motion is strongly forced* satellite (2) is "very stable" m roll-yaw 
but is "unstable" m pitch with forced motion; satellite (3) is 
"unstable" m yaw* "stable" m roll and pitch and pitch has forced 
motion* satellite (4) is similar to satellite (l) except that yaw and 
pitch are "destabilized" by the aerodynamic torque. 
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for satellite ( 2 ) with the 


( Compare , for example, the value of 
value for satellite (l) ) 

From the describing differential equations of controlled motion 
and the example trajectories of Figure 4 17 — Figure 4 20 it is apparent 
that for all initial state m S the station-keeping part of the con- 
troller can keep the state space trajectory from departing S by a signi- 
ficant amount (except, perhaps, when large unaccounted for disturbances 
overpower the station -keeping part of the controller). Since the minimum 
fuel expenditure of the station-keeping controller for one orbit is not 
known, the steady-state fuel cost per orbit obtained m the simulation 
runs cannot he compared with an absolute minimum However, since from 
Figure 4 17 — Figure 4 20 it is clear that the single-axis satellite 
motions are approximated, m at least a piecewise sense, hy the approxi- 
mation motions of Section A and this section, (Section B) the steady-state 
fuel cost is considered to be nearly a minimum Thus, the station-keeping 
part of the controller is considered to perform satisfactorily (in 
Chapter VI cost, error and other performance measures are evaluated for 
particular satellites of particular weights There it is found that the 
weight of fuel used m one year is very small compared to the weights of 
the satellites ) 
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V. DERIVATION OF ACQUISITION CONTROL LAWS 


A. EXTENSIONS OF BUSCH'S SOLUTION 

In Chapter III, Section D, Part 2 it was mentioned that Busch, using 
EMP and a reverse- time integration method, has found a nearly minimum- fuel 
optimal feedback control law which for a "stable" satellite results in 
acquisition to a region much larger (6.3 x 10 radians) than S (see e g. 
Figure (4.20) ) . The control law works well for x ^ 0.5 radians, 

1 = 1 , . . . , 6, but, for some large initial angles (about 6o°) , the coupling 
between the controlled roll motion and the controlled yaw motion via the 
controls through the pitch angle and the large value of cos ^0^ in the 
second of Equation (2.14) have a destabilizing effect. The Busch control 
law results m an acquisition time of about the time of one-half orbit from 
initial angles of about 25°. Acquisition from larger angles requires much 
longer acquisition times, and, m some cases when the satellite is ' un- 
stable" (e.g., satellite (3)) and the initial angles are large, acquisition 
takes a much longer time, (in some cases the motion can become uncontrol- 
lable so that acquisition cannot be accomplished. ) 

Thus, for the controllers of the satellites’ motions to perform 
satisfactorily, the Busch control law must be modified so that (l) the 
region S can be acquired even if large imperfections exist m the con- 
trollers, (2) the time of acquisition from x ^ 1 5 radians, l = 1, ...,6, 
is the time of one- quarter orbit or less, and, (3) "unstable" satellites 
such as satellite (3) are controllable for large initial angles. 

1. Via Phase Plane Techniques 

In component form the Busch control law can be written as 


= (-Nj/2)(SGN(x 2 |x 2 | + 2x J ) + SGN(x 2 + 0.1 x ± )) 


u r 



-N^SGN(x 4 ) 

0 


IF x^ > 1.7 X | x 3 | 
IF < 1.7 X |x 3 | 
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-J^SGN(xg) 



IP x| > 2 X jx 5 | 
IF Sg < 2 X ix,- | 


(5.1) 


where u_^, i = 1 , 2 , 3 ; is the nondimensional acquisition control torque 
measure numbers , I 1 ! = 1, and SGTT( ) is the signum function. Figure (5.l) 

shows sketches of the two types of control switching curves m the phase 
planes. It should he noticed that N 1 , l = 1, 2, 3, are all unity, and, 
therefore, the parabolic switching curves m the yaw and pitch phase 
planes are the same as for the minimum- fuel control of x" = u when 
I u y I = 1.0. Thus, if the magnitude of the control components are in- 
creased to shorten the acquisition time, an obvious modification to the 
Busch control law is the modification of the parabolic switching curves 
which makes them compatible with the larger values of IF , l = 1,2,3, i.e., 
m the equation x' = lax, the coefficient a (=211) is modified with IF 

Clearly, N = 1.0, l = 1,2,3, are too small to give acquisition 
from jx^l = 1.5 radians, l = 1,...,6, in a T-mterval of 1.57 i.e., in 
the time of one-quarter orbit. Consider acquisition from (x,x') = 

( 1.5 >1.5) for x" = u m minimum time It is not difficult to see that 

for |u | = 35 = 10 the acquisition time is ? - r = 0 95* From a 

nicix « x o 

study of the minimum- time acquisition of other simple systems i.e., 

2 2 2 

x" + a x = u, a > 0, a <0, it was found that minimum- time acquisition 
is accomplished m T f " t-q “ 1.57 if N was large enough. In particular 
the required values of IT ranged from a little larger than 2.0 to almost 
ten (depending on a ) . Thus, it seems that IF = 10, i = 1,2,3, are 
nearly lower hounds for the magnitudes of the control components for 
minimum- time acquisition. For minimum- fuel acquisition IF = 10, l = 1,2, 3* 
seem, certainly, to he lower bounds. Since values of IF , i = 1, 2, 3, 
which are much larger than ten (say, one- hundred) can result m very poor 
performance (such as higher cost and no acquisition to S ) when the usual 
imperfections m the controller (e.g. time delays) are present, the initial 

value assumed for each of N* , l = 1,2,3* was taken as ten. (in the next 

1 + ^ 
part it is seen that "spiraling m" to S from x ^0.1, i = 1, . .,6, 

should be avoided for low fuel cost. "Spiraling m" is avoided if 
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Figure 5.1 Switching Curves of the Busch Control law m the Taw and 


Pitch, Planes. 




Np i = l,2,3jr are large enough (^10 ), since the parabolic switching 
curves are nearly coincident with the "zeroing" part of a trajectory 
from a point on a parabolic switching curve with Ix^l 0.1, l = 2,k,6. ) 

Since Busch was successful in using a straight line for the other 

switching curve m the phase plane of the least "stable" motion and since 

the switching logic for such a curve is very simple, straight lines with 

various slopes were studied as possible switching curves. Phase plane 

methods with the constraints on the time of acquisition and fuel cost 

were used initially m the study m the hope of saving computer time. 

Since the initial values of the state variables considered were no larger 

than 1.50 radians, the control- on part of the phase plane trajectories 

2 

very nearly coincided with parts of parabolic curves such as x 2 = 2N^x^ 
if ^ 10. Thus, since the behavior of the phase plane trajectories 
was generally known when the trajectory was m a control- on region, the 
phase plane analysis was limited mainly to the control- off regions 

The phase plane analysis consisted of (l) the calculation of the 
slopes of control- off trajectories at numerous points m the phase planes 
by assuming the trajectories m the other two phase planes to be various 
nominal points and (2) for switching lines of various slopes the estima- 


tion of the time between control-on intervals with the aid of Ax = 

l 

(^i+l) AT off ’ 1 ~ 3>5, where S: is the (integral) mean rate value. 

Equations (2.lU) were used for |x | 1.0, l = 1, ...,6, and equations 


(2.17) and (2.18) were used for |x 


0.1, 1=1,.. ,6. The worst 


possible values S = jt/b and Q = ^/4(or -jp) were used m each case 
for each satellite. The results of this analysis can be summarized as 
follows (l) the shope of the switching line which gave a maximum control- 


off time interval of A-r __ 1.0 also resulted m most cases m the 

off 

lowest estimated fuel cost (the total control-on time for acquisition was 
estimated to be about 20-30% of the time of acquisition i.e., < ^' r on ~ 
0.3“0.5 *)j (2) if the motion was not very "stable" (e.g. satellite (l)), 
some of the smaller (m magnitude) slopes tried resulted m chatter motion 
which was very costly and time consuming (especially for |x_J < 0.01, 

1 = 1,3)5, where the sine forcing terms are most influential), and, (3) 
slopes of the switching line which where much greater (about lOx) resulted 
m a much higher cost estimate (about 200% higher). The slopes of the 
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switching lines which are ejected to result m a satisfactory acquisition 
control law are given m Table 5.1. These slopes are given m Table 5 1 
for each phase plane of each satellite and are generally compromises 
which are arrived at by placing the greatest emphasis on fuel economy and 
simplicity. For example, the shope of the straight line of the type 2 
switching curves is larger than was estimated as needed for large angles. 
The estimated values of the slope for pitch of satellites (l), (2) and 
(k) for large angles ranged from about -0.75 to -1.5, however, for small 
angles ( |x j < 0.01, l = 1,3,5) the magnitudes of the slopes were esti- 
mated to be about -6.0. Thifo , since pitch must be ” zeroed” faster than 
roll and yaw to avoid detrimental coupling and since other curves which 
give a variable slope are not as simple, a slope of -2 is a compromise. 

In the next part of this section the maximum principle is used 
to check and/or offer modifications to the control laws of this part. 

2. Via Pontryagm's Maximum Principle 

Using the reverse- time integration method, the maximum principle 
was applied to the acquisition problem. That is, the equations of motion 
(with the ‘coast function" optimal control) and the adjoint equations were 
written m backwards time and integrated with the aid of the digital com- 
puter. (For the program see Appendix F. ) The solutions of the equations 
of motion were plotted as phase plane trajectories (i.e., projections of 
the trajectories) so that the control switching points in the phase planes 
could be easily observed. 

Since S is small compared to the scale used m the phase plane 
plots and since, by "smoothing" the corners of S , the final (forward 
time) adjoint vector can span the six- dimensional vector space, the origin 

•J- 

of the state space was considered to be the final goal instead of S . 

That is, since the origin of the state space and S + are almost equiva- 
lent with respect to EMP and since computationally it is simpler to take 
x(T f ) = 0, the origin was acquired. 

The linearized equations of motion [Equations (2.17) with the 
aerodynamic torque terms included for satellite (4)] were used m this 
application of the maximum principle. The reasons for using the linearized 
equations are (l) the linear and nonlinear optimal solutions agreed 
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TABLE 5 1 MODIFIED VERSIONS OF BUSCH'S SWITCHING CURVES 
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very well for Jx^j < 0.1, 1 - 1,3,5, m the cases tested*, (2) for 
|x [ > 0.1, x = 1, 3,5, the nonlinear equations are highly coupled and 
result m erratic control switching points (since control switching for 
one axis depends on the motions of the other axes), and (3) the computer 
time for each solution run of the nonlinear system of equations was too 
great (about ten minutes) to be practicable (since many solutions were 
needed) . 

Some typical solutions of the linear optimal acquisition equations 
are given in graphical form m Figure 5.2. It should be noted that no 
trajectories "spiral m" toward the origin. "Spiraling m" is not a 
characteristic of an optimal trajectory m the case IT > 10, i = 1,2,3- 
Indeed, for such large values on W three switches m each control 
component is generally a maximum number m the time interval of " T 0 ~ 
1.57 since the adjoint variables on which the switches depend generally 
give (m backward time) a short control- on interval followed by a long 
control-off interval or with a short control-off interval and then a long 
control-on interval). 

In Table 5.2 are given some of the data obtained from the thirty- 
two backward time solutions. In particular Table 5.2 presents the range 
of values for the slopes of the switching lines which pass through the 
origin and the point at which the control switches from off to on. Also 
presented are the consensus values (for lines drawn through the greatest 
number of switching points) of the slopes - and the best values for stabili- 
zation of the motion. 

In the next section. Section B, the final states of the backward 
time runs are used as the initial states m the solutions of the full 
nonlinear equations of suboptimally controlled motion. These forward 
time solutions of the nonlinear equations of suboptimal controlled mo- 
tion and their fuel costs are compared with optimal linear solutions and 
their costs m the performance evaluation chapter. Chapter VI. 

B. ACQUISITION CONTROL LAWS WHICH PERFORM SATISFACTORILY 


*Also, Busch found very good agreement m many comparisons 
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TABLE 5 2 

SLOPES FOR THE STRAIGHT LINE SWITCHING 

CURVES AS 


OBTAINED FROM THE OPTIMAL (LINEAR) SOLUTIONS 

SATELLITE 

AXES 

RANGE OF SLOPES 

CONSENSUS 

BEST FOR 
STABILIZATION 

( 1 ) 

YAW 

-1 01 -> + oc* 

-1 00 

-1 01 

"STABLE" 



(OR -K) 4l) 



ROLL 

-1 28 -> + CO# 

+0 33 

-1 28 




(OR -1 13) 



PITCH 

-1 01 - -0 89 

-0 92 

-1 01 

( 2 ) 

YAW 

! -0 77- -0 59 

-0 69 

-0 77 

VERY "STABLE" 
IN ROLL AND 
YAW, "UN- 

ROLL 

-0 38 - -0 19 

-0 35 

-0 38 

STABLE" IN 
PITCH 

PITCH 

-2 1- -1 63 

-1 8l 

-2 1 

(3) 

YAW 

-0 86- + 

+0 42 

-0 86 

"UNSTABLE" 



(OR -0.67) 


IN YAW, 
"STABLE" IN 
ROLL AND 

ROLL 

- 0 . 91 — + co ~ 6 ' 

-0 61 
(OR 40.40) 

-0 91 

PITCH 

PITCH 

+1 03- +1 5 

+1 25 

+1 03 

00 

YAW 

_o 9 ^-— + ^ 

-0 77 

-0 94 

"STABLE" 



(OR 40 39) 


BUT GREATLY 
AFFECTED BY 
AERO 

ROLL 

-0 85 — + ccw 

+0 42 
(OR -0 84) 

-0 85 

■ TORQUE 

PITCH 

-1 33-* -0 89 

-1 02 

-1 33 


Depends on initial condition of adjoint variables , but, generally if 
yaw then roll —10 and vice versa 
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In Section A, Part 1 three types of phase plane switching curves 
(Table 5*l) were suggested, by the phase plane analysis as possible 
"improvements"* of the Busch switching curves In Part 2 of the last 
section the linear equations for an optimal acquisition solution were 
integrated m backward time. Some of the data optamed from these inte- 
grations are presented m Table 5 .2. 

If the consensus values for the slopes of the switching lines (Table 
5.2) are compared to the slopes given m the sketches m Table 5.1, good 
agreement is found m about half of the cases. The exceptions are pitch- 
satellite (l), roll- satellite (2) pitch- satellite (4-), which are each off 
by nearly a factor of two, and roll or yaw of satellites (l), (3) and (4), 
which are mostly off by a sign. If the best values for stabilization 
(the greatest slopes obtained from the optimal solutions) are compared 
to the slopes of Table 5.1, the agreement is found to be good m every 
case. 

Of the two types of switching curves, the parabolic curves have the 
smallest deviation from the optimal. This is, of course, expected for 
control magnitudes of ten or greater. From Figure 5.2 it is clear that 
the parabolic switching curves will work very well when }x |, 1 = 1, .,6, 

are small, however, for large angles and rates it is possible that some 
modifications will be needed. 

The following control law was used in the initial tests of acquisi- 
tion of the satellites from large angles [using the nonlinear equations 
of motion, (2.14) - (2.l6)] 

\ = (-r /2) {SGSr(x i+1 jx i+1 1 + 20x i ) + 

SGI r(x x+1 + QGm x - kj)]}, 1 = 1,2,3 (5.2) 

where M , the negative of the slopes, and the ” change- m- slope" factor 
Q^, 1 = 1,2,3, were as given m Table 5*3 (if Q as large, SGN 
(Q - |x |) is positive for all x and has essentially no effect ) 

Using some of the final conditions (backward time) of the optimal 

*Busch considered only one "stable" satellite. Here, consistent acquisi- 
tion to a smaller region from larger angles is obtained for four ("stable" 
and "unstable") satellites. 
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(linear) solutions for initial conditions, the full nonlinear equations 
with the control given by (5.2) and Table 5-3 were integrated with the aid 
of the digital computer. Since some chatter motion was expected along 
the straight switching lines m some planes, a time delay was built into 
the control. The time delay was mtially simulated by making the control 
a f unction of the value of the state at the end of the previous step of 
integration. The integration and the plotting (by the Stanford Cal- Comp 
Plotter) were carried out in two steps so that the small details of the 
phase plane plots could be observed. First, the digital computer inte- 
gration proceeded from |x.J « 1 5.> i = 1, ...,6, down to |x | ^ 0.02, 
i = 1,3,5. Then, using the final values of the first run, the second 

digital computer run was made down to jx [ ^ 2.0 X 10 \ 1=1,.. ,6. 

+ 

(Acquisition from S to S was simulated on the analog computer 
See Figures (4.17) - (4.20).) 

The initial runs (for the four satellites) were quite time consuming 
and only seven out of the twelve phase plane curves of each run appeared 
as expected. The trajectories m the phase planes of satellite (l) did 
not proceed toward the origin as directly as those trajectories of 
Figure (5.2). The reasons for this were considered to be (l) the slope 
of the switching line m pitch was too small (in magnitude) and caused 
the acquisition m the pitch plane to he too slow, and, (2) the slow 
acquisition m pitch resulted m harmful cross- coupling between the roll 
control and yaw and the yaw control and roll The fuel cost for this run 
was J = 15.87 and the time (t^, - t ) of acquisition was greater than 
2.0. The initial integration run of satellite ( 3 ) gave an unstable so- 
lution i.e., after a short time the state space trajectory began moving 
away from the origin and continued until there was no hope of acquisition. 
This unstable behavior proved to depend somewhat on the initial conditions, 
however, the main reasons for this bad performance are considered to be 
the same as the reasons for the poor performance m the initial run of 
satellite (l). The initial integration run of satellite (2) showed that 
early chatter motion m roll and yaw occurred. This early chatter motion 
appeared to be the reason for the cost and acquisition time being greater 
than expected. The (nondimensional) cost was J = 9*91 and the (nondi- 
mensional) time was 1.74. The initial solution for satellite (4) was 
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TABLE 5 3 THE PARAMETERS OF THE CONTROL LAW FOR THE INITIAL TEST 


SATELLITE 


( 1 ) 


( 2 ) 


( 3 ) 


00 


AXES 


YAW 


ROLL 


PITCH 


YAW 


ROLL 


PITCH 


YAW 


ROLL 


PITCH 


YAW 


ROLL 


M 

l 

Q x 

1.0 

10 0 

1.0 

10 0 

1 0 

10 0 

0.5 

10 0 

0 5 

10 0 

2 0 

10 0 

1 0 

10 0 

2 0 

o 05 

a 

0 05 

1 0 

10 0 

1.0 

10 0 




















































generally as expected except for a little early chatter motion m roll 
and yaw. The cost and time of acquisition compared very well with the 
optimal (linear) acquisition cost and time In all of the initial runs 
(except the run for satellite (3) which did not result m acquisition) 
the parabolic control switching curves gave good results with only a 
slight amount of chatter motion to S when m one case a "zeroing" 
part of a trajectory missed S + . 

The following changes were made m the phase plane switching lines 
of satellites (l) , (2) and (3), m the hope of improving the performance 
of their acquisition controls. (The initial run for satellite (4) indi- 
cated that its control resulted m satisfactory performance so that no 
changes m the control for satellite (4) were made.) The slope of the 
switching line for pitch- satellite (l) was decreased from -1.0 to -2 0 
m the hope of achieving quicker acquisition m pitch so that the detri- 
mental cross coupling between roll and yaw would be eliminated. The 
slopes of the switching lines m roll and yaw for satellite (2) were de- 
creased from -0.5 to -1.0 m an attempt to alleviate the early chatter 
motion m roll and yaw and, thus, reduce the fuel cost and time of acquisi- 
tion. The slope of the switching line m roll- satellite (3) was reversed 
m sign for (x^l ^0.05 (by changing 10 0) and was increased m 

slope from -2 0 to -1 0. This modification was made m an attempt to 
acquire [for the unstable satellite (3)] the region S (Since pitch 
for satellite (3) is so highly "stable", the negative slope of the switching 
line in pitch was retained. ) In the remaining simulation runs of the 
suboptimally controlled nonlinear system the computation time was shor- 
tened by removing the fixed time delay and introducing a more natural 
but slightly varying time delay. This ^was accomplished by introducing 
the stepsize cutting limiter of Part 1, Section A, Chapter IV. Por a 
stepsize of = 0 002, a time delay of about t^ = 0.125 sec was 
built into the control by limiting the number of cuts to four.* The ac- 
curacy was not significantly affected by this limit The step size was 
taken as Ar = 0.0002 for the second part of each run. This gave an ef- 
fective time delay of about t^ = 0 01 sec during this part The smaller 

*Recall that t, the nondimensional time (m radians) is equal to nt 
•where n is the average orbital rate and t is the real time m seconds 
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time delay was required when the trajectory encountered the region S + , 
since much larger time delays resulted m the overshooting of S + . The 
larger time delay was used during the first part of each run to save 
computer time. The large time delay was considered a worst case. 

The results of the second set of digital computer simulation runs 
for the four satellites with the new control laws [same as m Table 5-1 
except for roll- satellite (3)J and the new time delays showed definite 
improvements. Satellite (l) acquired the region S + within the time of 
one- quarter orbit and the fuel cost was reduced by about 25% . The re- 
sults for satellite (2) showed an improvement of about 10% in the fuel 
cost and a slight improvement m the acquisition time. Satellite (3) "was 
again -uncontrollable hut this time it took longer for the trajectory to 
begin moving away from S . New initial conditions (final conditions of 
an optimal (linear) solution) were used in the test of the acquisition 
control of satellite (4). Again, the results for satellite (k) were 
satisfactory. The fuel cost and time of acquisition compared very well 
with that of the optimal (linear) . 

The pitch control of satellite (3) was again modified. The positive 
slope of the switching line (when |x,_ [ > 0.05) was changed to a negative 
slope m the hope of reducing the large angle coupling between the yaw 
control and roll and the roll control and yaw by rapidly reducing the 
pitch angle. Simulation of the acquisition of satellite (3) with the 
newly modified control was made for several initial conditions. These 
simulation runs gave satisfactory results although the fuel costs were 
35 - 67 % greater than the optimal (linear) fuel costs. 

In Figure (5-3) - Figure (5.6) are shown phase plane plots of the 
phase plane projections of one of the worst case acquisition trajectories 
for each satellite. These are worst case trajectories since they show 
more chatter than others and since the fuel costs and/or times of acquisi- 
tion are generally greater. The right plot on each page of Figures (5*3) - 
(5*6) is an enlargement of the origin. Of course, these plots do not 

*"h r 

show all of the acquisition. Once a phase plane projection of S l see 
Figure (k 17) - Figure (k 20)] is acquired by the phase plane projection 
of the acquisition trajectory, the station- keeping part of the controller 
[see Figure (6. l) ] takes over and performs the acquisition form near the 
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Large Ysi* Angles 1 Small Yaw Angles 

Figure 5 5 Suboptimally Controlled Acquisition Motion of Satellite (3)^ 

J - 9*75; fjj, - t = 1.63, 9_ = it/4 
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projected boundary of S + to the projection of S. (Since two widely 
different levels of control torque are needed for satisfactory performance 
m acquisition and station-keeping, the region S was acquired by the 
acquisition part of the controller. Otherwise, rapid chatter motion and 
wasted control effort result from the control torque being too large for 
the region of station-keeping. See the data of Chapter ¥1 on the fuel 
cost versus the size of N , i = 1,2,3, and the size of S.) 
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VI. COMPLETE ATTITUDE CONTROL AND PERFORMANCE EVALUATION 


[The complete attitude control system consisting of the station- 
keeping controls derived m Chapter IV and the acquisition controls 
derived m Chapter V are considered to give satisfactory performance 
A block diagram of this complete attitude control system is given in 
Figure (6.1). 

In the remainder of this chapter the performances of the two parts 
of the complete control system are evaluated m terms of the performances 
of other systems and m terms of the weight of fuel required for a year 
of control as compared to the weight of the satellite. Also, presented 
here are performance limitations imposed by imperfections. 

A. ACQUISITION CONTROL COMPARED TO OPTIMAL (LINEAR) ACQUISITION CONTROL 

In Section B of Chapter V, the selected acquisition control system 
was tested and found to perform satisfactorily, however, a thorough 
comparison of the fuel costs found in the test runs with the optimal 
(linear) fuel costs was not made. Although the optimal (linear) fuel 
costs are probably very conservative, they are considered to he a lower 
hound which is large enough for the comparison to he meaningful 

The average nondimensional fuel cost for the acquisition simulation 
runs of satellite (l) was J = 12.4. For satellites (2), ( 3 ) and (4) the 
average fuel costs were 10.2, 14.1 and 13*2, respectively. These average 
costs differed from the average optimal (linear) fuel cost by 34%, 15 %, 
43% and 36 % for satellites (l), (2), ( 3 ) and (4), respectively. The 
percent differences of the average costs seem quite high, however, con- 
sidering the strong cross coupling which occurs when the angles are large 
( ^ 6o°) , the fact that a simple fixed switching logic control is used and 
the fact that the optimal (linear) cost is not a greatest lower hound, 
these percent differences are not unsatisfactory. 

Several simulation runs of the nonlinear acquisition system were 
made for small angles (~ 5°) from the final conditions (backward time) of 
optimal (linear) runs. The percent differences for the small cases were 
about 15 % . 

The variation m the percent difference for each individual run was 
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NOTES : (l) The index 1 = 1,2,3 for yaw, roll and pitch, respectively 

(2) = 1.0, Mg^lO and !^ = 20 

( 3 ) The distance d (see Section B, Chapter 4) depends on the 
satellite, the strength of station-keeping control and on 
the controllers time delay 

(1) U , 1 = 1,2,3, are the nondimensional control components 
(either acquisition or steady-state) 

Figure 6 1 Block Diagram of the Complete Attitude Control 
System 
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from -4.1 % to +62%. The extreme percent differences belonged to the 
large angle runs. The suboptimal system had an advantage of 4 1% over 
the optimal system m one of the runs for satellite (4), however, the time 
of acquisition for the suboptimal was almost 20% greater than that of 
the optimal. The maximum percent ldfference of 62% was for the run of 
satellite ( 3 ) which was given m Figure (5.5). 

One of the mam reasons for the greater fuel cost of the suboptimal 
system is the chatter motion whri ch occurs during acquisition from some 
initial conditions. This chatter motion can be avoided, even with a fixed 
switching logic control, by increasing the values of M (the negative 
values of the slopes of the straight line switching curves ) . If the 
values of are increased until no (or little) chatter motion occurs, 

the time of acquisition generally decreases, but, the fuel cost increases. 
However, unless the slopes are increased until no (or little) chatter mo- 
tion occurs (about a factor of ten), the costs do not generally increase 
by more than 100% . 

The nondimensional fuel cost values can be translated for particular 
satellites into the weight of fuel used m, say, pounds. If, while a gas 
jet is on, it is assumed that the gas flow rate and exhaust velocity are 
constant, the weight- flow for a gas jet couple is given by 

r I n 2 

w = ~ 2l 2T~ I ' 1 = (6-D 

1 

where is the moment arm distance, g is the acceleration of gravity 

and x> is the magnitude of the exhaust velocity. If both sides of (6.1) 
are multiplied by n ^Ar (recall that t = nt), the result is 


w = w At = ST'Ar 
11 1 



( 6 . 2 ) 


which is the weight of fuel used in the time interval, At, that the 
gas jet couple is on. 

The ratio u/g, whose reciprocal is a factor m (6,2), is usually 
referred to as the specific impulse, I . A realistic value for v is 
1500 ft. /sec. so that a reasonable specific impulse is I = 46.5 sec. 
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— p — *1 

In this case g/'O = 2.15 X 10 sec. . The values of n for the satel- 

lxtes considered m this investigation range from 0.95 X 10 to 1.13 X 

-3 

10 . In the following calculations of the weight of fuel used, the 

worst case value of n = 1.13 X 10 is used. The typical satellites 
referred to are those of Section D, Chapter II. Since the typical 
satellites are sperical or cylindrical and since the lever arm distance, 

2 is usually half the length of an axis, the ratio 1^/2^^ can he 
approximated by l/5 mf^, if the satellite is spherical, by l/4- for 
the longitudinal axis and by l/6 mg for the transverse axis, if the 
satellite is cylindrical, (m is the satellite's mass. For satellite 

(1) m = 1,860 slugs. For satellites (2), (3) and (4) m is 1550 slugs, 

1550 slugs and 310 slugs, respectively.) 

Now the nondimensional average acquisition fuel cost for each satel- 
lite can be translated into pounds weight. The values of J (given 
earlier m this section) for satellites (l), (2), (3) and (4) are 12.4, 

10.2, 14.1 and 13-2, respectively. Therefore, the approximate (worst) 
weights of fuel used on the average for acquisition by satellites (l), 

(2) , (3) and (4) are 0.88 lbs. (f x = 7 ft.), 1.28 lbs. (£ = 15 ft.), 

1.78 lbs. (i^ = 15 ft.) and 0.10 lbs. (i =45 ft.), respectively 

B. STATION-KEEPING CONTROL 

In Part 4, Section B of Chapter IV the station-keeping control sys- 
tem was simulated on the analog computer and found to perform satisfactorily. 
In this section the performance of the station- keeping control is eval- 
uated by comparing it to the hybrid station- keeping control system of 
Busch (whi'di contains a reaction wheel) and by evaluating fuel costs 
changes with various parameters of the system. 

1. Comparison with Busch's Solution 

Busch obtained station- keeping control m roll and yaw simply 
by leaving small (6.3 X lO - radians) circular regions of no control 
about the origin in the roll and yaw phase planes of the acquisition 
phase plane switching logic. This performed well for the "stable” 
satellite considered by Busch. (The steady- state cost per orbit for 
roll and yaw was about J = lCf .) 
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Since the pitch motion was forced by a sinusoidally varying 
term and since neither of the switching curves were of a "stabilizing" 
type, Busch's pitch acquisition control could only achieve complete 
acquisition m pitch periodically and the error m pitch grew periodically 
to about 0.1 radians. The fuel cost for this periodic reacquisition m 
pitch was high. To gam greater accuracy m pitch and to reduce the 
weight of the fuel used, Busch supplemented the pitch part of the control 
system with a reaction wheel. 

If the environmental torques acting on the satellite are nearly 
sinusoidal, the power required for "pushing" against the reaction wheel 
is nearly zero, so that, except for the weight of the reaction wheel and 
supporting components the weight of the pitch part of the controller is 
negligible. However, since the environmental torques are not sinusoidal 
hut have the character of a nearly constant forcing term with some periodic 
variations, the reaction wheel must be continuously "pushed" against m 
such a way that the wheel speed is continously increased. Since Busch used 
gas jets to produce the torque to hold the satellite steady while slowing 
the reaction wheel once it reached its saturation speed, his steady- 
state fuel consumption was increased by 21.9%* (This increase was due 
to slowing the reaction wheel when only cross coupling caused the wheel 
speed saturation i.e , all forcing terms were considered as purely sin- 
usoidal. ) 

The weight of the fuel used for station-keeping can he compared 

(for the same size satellite) to the weight required by Busch's system 

for station- keeping control. Consider the "stable" satellite (2) Busch 

-5 

suggests a reaction wheel moment of inertia of I =10 X I = 1.21 
2 w 3 

slugs- ft. . If the reaction wheel is a brass cylinder with a diameter 

of 1.5 ft. and a height of one foot, the weight of the reaction wheel is 

138 lbs. From Table 4.1 it is found that the cost of the steady- state 

-2 

pitch control for satellite ( 2 ) is J = 6.03 X 10 per orbit. The 

o Bj 

cost for a year (about 5-6 X Hr orbits) of pitch steady state control 
is J = 337* From Equation (6.2) the weight of fuel used m pitch for 
a year of station- keeping is found to he 8.47 lbs. (i = 3 ft.). In the 
same manner the weight of the fuel used for roll-yaw station- keeping is 
found to he 1 13 lbs. (i = 15 ft.) per year which is about the same as 
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Busch's roll- yaw control used for the 50,000 Tb. satellite. 

Clearly, the weight of the fuel used for station-keeping by the 
control of Section B, Chapter IV is much less (about 9^) than the weight 
required by Busch's station- keeping control system. 

2. Fuel Expenditure Per Orbit For Values of Various Parameters 

The fuel cost for station- keeping varies with the values of 
parameters such as eccentricity (e), inertia parameters (k ), strength 
of the control (IT^), size of and the peak voltage of sensor 

noise. Simulation runs (as m Section B, Chapter IV) were made to de- 
termine the effect of the variation of these parameters on the fuel cost. 
(The results of the runs also offer a check on the behavior predicted 
for the satellites by the analysis of the simple motions of Chapter IV.) 

To test the effect of eccentricity on the fuel cost, simulation 
runs were made for the inertia parameters of satellite (3) which is 
"unstable" m yaw but "stable" m roll and pitch. From these runs it 
was found that the roll- yaw fuel cost varied by only 3% for a variation 
in e from 0.0 to 0 05, however, the pitch fuel cost varied directly as 
the eccentricity Figure (6.2) shows the pitch fuel cost as a function 
of eccentricity. 

To test the effect of the inertia parameters on the fuel cost 
per orbit, five sets of values of k^, i = 1,2,3, (which correspond to 
five satellites whose earth- pointing motions range from highly "unstable" 
to highly " stable" ) were used m the simulation runs . The change m the 
pitch fuel cost was negligible when e = 0.05 was used. A second set of 
runs was made with e = 0.01 m the hope that the reduction of the size 
of the relatively large forcing term would result m noticeable changes 
m the pitch fuel cost with changes m k^, however, the results were 
the same for this set of runs. The changes m the roll- yaw fuel cost with 
changes m and k^ were quite significant (a maximum increase of 

172%) These changes are represented m the chart of Figure (6.3) for 
®1 = E 2 = 0.01, e = 0.05 and d = 3 X 10 -5 . 

Satellites (2) and (3) were simulated m the investigation of 
the effect of the strength of control on the fuel cost. (These two 
satellites were chosen for this investigation since satellite (2) is 
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Figure 6.2 Pitch Fuel Cost Per Orbit vs. Eccentricity, N = 0 035, k = O. 99 , d = 0 
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of k^ and k , 


’’■unstable" m pitch and satellite (3) is "unstable" m yaw.) The dis- 

-5 -5 

tance d was taken as 2.0 X 10 ; 2.5 X 10 m the pitch control and as 
-5 -5 

2 X 10 and 4.5 X 10 m both the roll and yaw controls. Recall from 
Chapter XV that m pitch, when the sinusoidal forcing term is most sig- 
nificant, the distance d should be large enough to insure a control- 
off interval m any interval which contains two consecutive control 

swatches to iEL. Also, recall that for "stable" single-axis motions 
o 2 

(e.g those of x" + a x = v) the fuel cost should decrease with increases 

m a(-10~ . For "unstable" single-axis motions (e.g. those of x" + 

2 2 „ 1 |_ 
a x = v, a < 0) the fuel cost should increase with d(^10~ ), however, 

d must be just large enough to insure a control-off interval in any 

interval which contains two consecutive control switches to i'N 

Figure (6.4a) shows the roll-yaw fuel cost per orbit of satellites 

(2) and (3) plotted as a function of the control strength for the various 

values of d. Since the pitch fuel cost varies insignificantly with the 

inertia parameters, only the pitch fuel cost per orbit of satellite ( 3 ) 

is shown plotted m Figure 6.4b. The points on the curves of Figure 6.4b 

marked "overshoot" correspond to the values of SL (for a time delay of 

about 0.5 sec.) for which control- on pitch trajectories overshoot the 

small strips (of width d) m the pitch phase plane. For 1'!^ greater 

than these values the values of d are not large enough to insure a 

control-off interval m any interval which contains two consecutive 

control switches to idy 

The effect of the size of the region S on the fuel cost per 
orbit was investigated for satellite ( 2 ). The results of this investiga- 
tion are given if Figure 6 . 5 . The point (0.0,0.0) aided m the con- 
struction of the curve for roll- yaw. (This was possible since the origin 
of the roll- yaw state space is an equilibrium point.) 

To test the effect of sensor noise on the error and fuel cost 
per orbit a low frequency (100 cps) Gaussian noise generation was used. 

The noise was added to the state variables as they entered the controller. 
Runs were made for satellite (3) with I'X L = 0.04 and d = 2 5X 10”\ 

The peak noise voltages used were 0.05, 0 1 and 0.2 times the peak state 
variable voltage of 1.0 volt. The error and fuel cost of each run with 
noise were compared to the error and fuel cost of the noiseless run. 
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Figure 6.1 Fuel Cost Per Orbit vs Control Strength. 
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The fuel cost was found to increase slightly (less than 2% ) while the 
error was at most 0.5 X 10 ^ radians or 5%. 

C. IMPERFECTIONS IN THE MODEL AND THE ULTIMATE ABILITY TO EARTH- POINT 

1. Nonrigidity of the Satellites 

The mathematical satellite models were made on the assumption 
that the satellites are rigid. This, of course, is not precisely the 
case, however, for the satellite configurations considered m this in- 
vestigation, the effects of the various causes of nonrigidity can he 
made acceptably small Contributors to nonrigidity include the slight 
deformation of the satellite structure m the small force environment 
(The largest gas jet required for acquisition gives only 0.17 Its. of 
thrust.), the motion of the valves of the gas jets and the gaseous fuel 

In the steady- state mode, the deformation of the satellites 

should be insignificant when, for example, an attached camera is to be 

-A 

aligned very accurately for earth- pointing to within 10 radians. For 
the satellites considered here the deformation is generally much smaller 
than the best machining tolerance available today and is, therefore, 
insignificant 

The motion of the gas jet valves is very rapid so that the 
valves are opened or closed m just a few milliseconds. Thus, the gas 
valve must attain a speed of the order of 100 m./see from rest (with 
respect to the satellite). After attaining this speed the valve must 
again be brought to rest;, The forces required to accelerate the values 
and to bring them to rest can cause a significant torque (about 0.1 ft.- 
lbs.) unless care is taken m their design. 

Since the pressure due to a completely gaseous fuel is (for all 
practical purposes) uniformly distributed over the inner surface of the 
fuel container, the torque due to the pressure forces is essentially zero 
(regardless on the containers shape or its position m the satellite). 
Therefore, the third contributor to the satellite's state of nonrigidity 
has an insignificant effect on the satellite's earth-point mg performance. 

2. Gas Jet Misalignment 

Each gas jet pair was assumed to be aligned such that the torque 
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of the couple was "about" a principal axis of inertia for the mass center. 
In practice the orientations of the principal axes are known to within 
only about 0.01 radians. Hence, when a gas jet pair is on, it is probable 
that a torque with a magnitude of about 0.01 of the gas jet pair's torque 
is applied about the other principal axes. It is most likely that these 
misalignment torques will not improve the control systems performance, 
however, their effects on the accuracy are considered to be insignificant 
but they can result m as much as a 12% increase m the fuel cost. 

3. Sensors and Error 

In this investigation it is assumed that the state of the at- 
titude motion is available from sensors. In the last section. Section B, 
it was found that low frequency (100 cps) noise added to the s bate 
variable signals from the sensors has almost no effect on the performance 
(even if its peak voltage is 0.2 of the sensor's output voltage). How- 
ever, if there is a bias error of 10%, say from the misalignment of the 
sensors, the earth-pointing accuracy will suffer by 10% at times and the 
fuel cost will generally increase. 

4. Others 

The mathematical satellite models used m this investigation do 
not account for very large forces such as those due to the motion of a 
man on board and the collision of the satellite with a very high momentum, 
but "non- fatal" meteoroid. 

Calculations based on Ax =* i = 1,2,3, and on 

n I Urn Ar 55 i! m Ap i = 1,2,3, where Nm is the effective value of the 
ith component of the nondimensional meteoroid torque, m is the meteor- 
oid's mass, is the moment am for the i^ 1 axis and Ap is the change 
m the magnitude of the meteoroid's velocity, showed that (for Ap be- 
tween 20,000 ft. /sec. and 60,000 ft. /sec. and for i between 1.0 ft. 

and 15 ft.) the station- keeping control should he able to easily accom- 

-5 

modai e meteoroids with masses up to about 5 X 10 lbs. several times per 
orbit. (Accommodation of meteoroids of this size is considered to occur 
rarely. See Section E, Chapter II.) The acquisition control, which acts 
like a hack-up station- keeping control with some decrease m accuracy, 
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should be able to accommodate meteoroids up to about 0.05 lbs. (with 
little loss of accuracy) . 

A man m motion aboard the satellite can exert forces on the 
satellite which cannot be compensated for by the station- keeping control. 
(The force on the satellite due simply to the man casually raising his 
arm above his head is about 10 lbs.) The acquisition control cannot 
compensate for these forces unless (for the satellites considered m 
this investigation) the nondimensional control magnitudes, N , 1 = 1,2, 3 > 
are increased to about 10,000 and the time delays are reduced to about 
10 microseconds. (Of course, if the satellites are much more massive 
than those considered here, say about 1000 tons, the values of H 
of the time delays do not need to be so extreme.) 
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VII. CONCLUSION 


A feedback control system for efficiently controlling the attitude 
motions of satellites m elliptic orbits about an oblate earth with an 
atmosphere was devised The criteria used for efficient (or "satisfactory") 
performance were (l) high- accuracy (10~ ^ radians) m earth-pointing 
after the acquisition of the earth- pointing mode has been accomplished 
within the time of one-quarter orbit from large angles (>6 q°), (2) 
minimum fuel expenditure, and, (3) practicality of the system. Although 
four particular satellite configurations ("stable" and "unstable") were 
assumed m the derivation, the devised control system performed well for 
a wide variety of satellite shapes, orbits and parameters of the control. 
(See Chapter VI.) 

The derivation of the control system proceeded m two steps. First, 
in Chapter IV the station- keeping part of the controller was devised. 
Pontryagm’s maximum principle, the "jump conditions" and the guidelines 
obtained from the minimum- fuel station-keeping controls devised for 
single-axis systems were used in this derivation. The maximum principle 
was applied to (l) the minimum- fuel problem which was considered as a 
problem m the theory of optimal processes with bounded phase coordinates 
and (2) the minimum- fuel problem with an integral constraint on the state 
(of the attitude) for maintaining high-acenracy earth-point mg. 

In Chapter V the acquisition part of the controller was devised. 
Pontryagm ' s maximum principal and phase plane methods were used to extend 
the Busch acquisition control to give (l) acquistion for a wider range of 
satellites (including "unstable" satellites), ( 2 ) higher accuracy and ( 3 ) 
acqusntion from larger angles an the more realistic time of one- quarter 
orbit o 

Two avenues of future research related to this investigation are 
considered to be (l) the theory (sufficiency conditions) for the long- 
time optimal control of highly forced systems with bounded phase coordinates 
and (2) the high-accuracy and efficient control of the attitude motion of 
a manned satellite (by, perhaps, devices other than gas jets). 
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APPENDIX A. THE GRAVITATIONAL TORQUE 


DeBra [ll] has shown that the gravitational torque for B* (the c m 
of the satellite) due to an oblate earth is suitably given by 



t (5qJ r E 2 /r 5 ) n^ x I R (l - 7 cos 2 y) 

' (^Jh//r 5 ) H x I N 

+ (!OqJr E 2 /r 5 ) cosy (n^ (A.l) 


where J = 3/2 (i^ - Igj/m^rg 2 = 1 63 x 10 q = nigG, I and I R are the 
polar and equatorial principal moment of inertia for E^ (the c m of the 
earth) nig is the earth's mass, r p is the earth's equatorial radius, G 
is the universal gravitation constant, np-j. is a unit vector directed 
from E* to B*, I is a unit vector which is parallel to the earth' s axis, 
7 is the angle between N and hp-,, and Ipq and Ip are second moment vec- 
tors of the satellite, B, relative to B* for nj^ and N, respectively 
Kane [22] shows that I R , for example, can be written as 



■l 



I. n 
ij -3 


(A 2) 


where n^, 1 = 1 , 2 , 3 , are mutually perpendicular unit vectors, b = N . 

n , i=l, 2 , 3 , and I , 1 , j = 1 , 2 , 3 , are moments and products of 

^ 1 J 

inertia of B relative to B* for n and n 

-1 —J 

Suppose n , 1 = 1, 2, 3 are parallel to principal axes of inertia 

of B for B* (see Figure 2 l) Then I =0 for 1 4 j. If n , 1 = 1, 2, 

1 J 3_ 

3, are right handed unit vectors and a^ = ng 1 • n^, 1 = 1 , 2 , 3 , then 
with the aid of A. 2 and a similar relation for Ip^ the terms which ap- 
pear as vector products m A.l can be written as 
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Sr, x 5r , * a 3 a 2 (I 3 - V Si + a l a 3 (I 1 ' 1 3 ) s 2 


+ (r 2 - - n 3 


I X I N = b 3 b 2 (I 3 - l 2 ) n x + b 1 b 3 (I x - I 3 ) n 2 
+ b 2 b l " I i^ -3 

—R 1 X — R ^ a 2 b 3 I 3 " a 3 b 2 I 2' > -1 + ^l 1 ! ' -2 

+ (a i b 2 I 2 “ a 2 Vl ) —3 

N Xl R = (b 2 a 3 I 3 - 'bgagig) 2 ^ + - t^a^) £ 2 

+ ^Ws " b 2 a l I l ) —3 (A 3) 


where 'I -“'I , i = j = 1, 2 , 3 

l ig 

Row the gravitational torque as given by A 1 can be written with 
the aid of equations A. 3 m the desired component form A more useful 
expression for the torque is obtained if a , and b^, l = 1, 2 , 3 ) are 
written as functions of the attitude angles 0 , l = 1, 2 , 3 ? and of the 
parameters of the orbit of the satellite, S, 9 , 9 (see Figure A.l) 

This will be done before equations A. 3 are substituted into A.l. 

It is not difficult to conclude with the aid of Figure A.l and 
trigonometry that 

R = smB [ sm (9 + 0^) n^ + cos (0 + 0^) n^ ] + cosSn^ (A 4) 


where 


is the unit vector given above, np^ 


is a unit vector 


which is normal to the orbit plane with its sense given by the right 
hand rule for increasing 6 , is a unit vector which together with 

and np> 3 form a set of right-handed mutually perpendicular unit 
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Figure A 1 Unit Vectors and Parameters of the Orbit 
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vectors and 0^ is the angle between the ascending line of nodes and 
the line "between E* and the perigee point of the orbit. 

Nov, if ng , i = l, 2, 3, are vntten m terms of n^, i = 1, 2, 
3, and if these expressions are substituted into A A, then the coeffi- 
cients of n , l = 1, 2, 3> in the resulting expression vill be b_^, 
1=1, 2, 3, respectively The coefficients of n , l = 1, 2, 3) m 
the expression for rijj are a^, i = 1, 2, 3> respectively 

If the three-axes Euler angles, 0^, 0g, 6^, are used as m 
Busch [6], it can be concluded that 


\ - °2°3 Si '°2 S 3 
\ ’ (s l s 2°3 + V3 5 % + ( °1°3 ' 6 1 S 2 S 3 ) 
£e 3 “ (s l s 3 - -1 + (S 1°3 + C 1 S 2 S 3 ) 


E s +s aS 3 


-2 ” S 1 C 2^3 


—2 + C 1°S*J 


(A 5) 


where c = cos 0 , s = sin 0 , l = 1, 2, 3 
i i i i 

From the first of A 5 it is observed that 


a l = 




a 2 = - C 2 S 3^ 


a 3 


(A 6) 


and from A h and A 5 it is observed that 


h - Ve° 2 c 3 + V 0 V 3 ' c 6 c i s 2°3 + Ve s i s 2°3 + °s s i s 3 


b 2 - VeV3 + Vi°3 • s B s e° 2 s 3 + C B° 1 S 2 S 3 • Ve s i s 2 s 3 
b 3 * Vl°2 + S B S 8 S 2 - S 8°S S 1 C 2 


(A 7) 


where, for example, e c = cos 6 and s_. = sm (0 + 0 ) 

O C7 p 

Nov with A. 3} A 6 and A 7 the gravitational torque, which is given 
by A.l, can be expressed m a form which is convenient for studying its 
effect on the attitude motion of earth- pointing satellites The sub- 
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stituuion of A. 3 into A, 7 results m 
where 


T = 1 n.. + T _ n_ - 1 - T _ n_ 

g gl -1 g2 -2 g3 -3 


T gl = (*3 " I 2 ) (3a 3 a 2 + 5 j [ t ^/ t ] 2 [1-7 cos 2 y3 

- 2J [r^r] 2 + 10J [r^/r] 2 cosy [a 2 b 3 + a b g ] } 

T g2 = (}J./r 3 ) Cl-j_ - I^) (3a 1 a 3 + 5J [r^r] 2 [1-7 eos 2 y] b .^ 3 . 

- 2J [r^r] 2 b 1 b 3 + 10J [r^/r] 2 cosy Ea^ + a^] } 

T g3 - (ii/r 3 ) (I 2 - I x ) [3a 2 a 1 + 5J [rjr] 2 [1-7 cos 2 y] 

- 2J [r^/r] 2 b^ + 10J [rg/r] 2 cosy [a^bg + a^] } (A. 8) 


In the terms of equations A. 8 the factors which involve only 
and b , l = 1., 2j 3 ? are given by 


a 3 a 2 — ^”2^2^3 ? a i a 3 ~ ^gSgC^ a 2 a i ~ “ e 2 ^2^3 


Vs * Ws'lSh + S 5 2 °0 3 0 C X S 2 O 3 ' S S 2 °0°lV2 C 3 

. 2 „ 2 / , 2 lu 

+ C S C 1 1 C 2 C 3 " V& 0 1 2 S 3 * ( terms in s x > * * • » \ ) 


V3 = VsVl^S + *bW B Z°3 - S 6% C 0 S 1 C 2 2c 3 


, 2 2 2 .. 2 3\ 
S & 5 C 0 C 1 C 2 S 3 “ e S C 1 C 2 S 2 C 3 ' terni3 an \ t s x ) 
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b 2^1 s 5 S e C 0 C l C 2 C 3 + S 5 C 5 S e S l C 2 C 3 ~ S 5^ S 0 C 2 3 e 3 S 3 

2 2 2 2 2 2 S 

+ s S C Q c l C 3 S 3 ~ VsVl S 2 C 3 + ( terms xn S 1 » **•> \ ) 

2 1l 

a 2 b 3 + a 3 b 2 = C 5 e i C 2 ( S 2 ” C 2 S 3^ + ( terms in ^ ) •••> ^ ) 


a 3 b l + ^ + 2 - Wl^S 


2 ^ 

+ (terms xn s 3 s J ) 


a l b 2 + V 5 ! - WlV^ • 2 Ve 0 2S 6 3 + °6 S 1°2 C 3 2 


2 *3 

+ (terms xn , s^ 0 ) 


From A =4 xt xs found that cosy = N 
cosy = s 5 s 0 


xip xs gxven by 


(A. 9) 


(A. 10) 
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APPEMDIX B. THE AERODYNAMIC TORQUE 


If a satellite is m orbit 100 miles or more above the earth's 

surface, then the ratio of the speed of B' x ' to the most probable random 

3 

speed of the air molecules is of the order of 10 . At these altitudes 
the mean free path of the molecules greatly exceed the maximum linear 
dimension of the satellite Hayes and Probstem [l8] have shown for 
bodies m such free molecular flow that the pressure and shear is given 
by 

P = (2 - f^pV 2 cos ty 
T = f^pV 2 Sin <1* cos 4 1 


f^ are the normal and tangential accommodation coeffi- 

_ | v P/ a | W here v^ 3 


is the 


where f and 
n 

cients, p is the density of the air, V 
velocity of a point P of the element of surface area of B relative to 
the free -stream velocity of the air near the satellite and i{j is the 
angle between the vector f>/ a and a line normal to the surface element. 
For the satellite attitude motions considered m this thesis the 


angular velocity of B m inertial space is such that error m taking 
v^ 3 = v^ is at most 0.0001$. In the following derivation V will 


be assumed to be 


B*/ s 


and 


The accommodation coefficients are defined by f. = (t - t )/t 

t i r i 

f = (p - p„)/(p - p. ) where, for example, p. is the normal 

XX Jl 3? 1 D 0 


momentum component of the molecules which are re-emitted from tne surface 
with a Maxwellian distribution at the surface temperature and the 

subscripts i and r denote incident and reflected. Hayes and Probstem 
conclude from the results obtained by experimentalists that f lies 
m the range between 0 8 and 1.0 and f is about unity 

If n 

— n 

area, say ds, and is directed out of B and if 

B*/a „ 

v ' = Vn , 


denotes a unit vector normal to the element of surface 

n is such that 

then the force on ds is given by 


dF 


= -pV 2 cos i|i[(2-f )cos n - f,. sm n, ]ds 

II “XI “*C 


(B 1) 
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where is a unit vector m the direction of the projection of -n 

onto ds . 

For a spherical satellite the integration of (B.l) over the half 
of the surface m the flow gives 

F - -CpSpV^n (B.2) 

where CL, - 2 - f - f . and S is the projected area so that the force 
D n t B w 

vector is parallel to v ' and is a drag force. 

If f = = 1 m (B.l); then c!F is parallel to v for an 

element of a satellite of arbitrary shape. Thus, the sum F will he a 

drag force only for any satellite if f = f ' = 1. 

n "u 

If f = 1 and f, = 0.80; the lowest value given above; then m 
n x 

(B.l); the term m parentheses can be written as 
n + 0.2 sm i|»(cos - sm i|rn) 

where n^ is perpendicular to n and is directed maximally away from 
ds'. Examination of this term shows that a "lift" component of the force 
is possible for some orientations of some satellites. For a cylindrical 

satellite with either circular or elliptic cross section m a nearly 

-3 

earth -pointing orientation; say 10 radians ; the term which contributes 

_3 

to the lift is smaller than 10 times the drag component for both 
lateral and end surfaces so that the "lift" component of the resultant 
force will be insignificantly small. For all other orientations of any 
satellite the ratio of "lift" to "drag" is still less than 0.12. 

The aerodynamic force used m the calculation of the aerodynamic 
torque for B* will be that given by (B.2). 

Kmg-Hele [23] considers the error to be less than 5$ when is 

taken to be 2.2 for both spheres and cylinders (i/d > l) if is 

based on the mean area; S; perpendicular to the direction of motion. 

The atmospheric density p has been determined experimentally and 
theoretically m resent years for altitudes above 100 miles. It has 
been found (see King-Hele [23] and Johnson [21]). that p decreases 
almost logarithmically with altitude between about 150 miles and about 
500 miles. If this be the case then 
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p = p exp[-K(r-r )] (B 3) 

Jr Jr 

where the subscript p denotes perigee, K is a constant to be deter- 
mined from the data and r is the distance from E* to B^„ The 
atmospheric density data is usually plotted on semi-log graph paper so 
that K can be read from the plots directly From these plots it is 

observed that K should be about 0 02 miles “ from an altitude of 200 

-1 

miles and should be about 0 01 miles from an altitude of 500 miles 

Also, it is observed that due to changes m the sun’s influence between 

day-time and night-time, the density changes by a factor of about two 

at an altitude of 200 miles and by a factor of about ten at an altitude 

of 500 miles. The maximum values occur m the daytime. If it is 

desired to fix p and K for any one orbit without, regard to day and 
Jr 

night, then for the most nearly correct values of p over the range of 
altitudes of the orbit and over the time interval of days it is best to 
pick the day-time value of p^ and the value of K which corresponds 
to day-time values of p near perigee. This is done m the thesis for 
a height of perigee of 200 miles In this case the density, p, as 
given by (B 3) is considered to be too large during night-time by a 
factor of two at 200 miles and by a factor of ten at 500 miles. It is 
considered to be too small by a factor of 0,7 during day-time at ^00 
miles and to be correct during day-time at 200 miles 

The torque of the aerodynamic forces acting on the elements of 
surface area of a spherical or high-accuracy earth-pointing cylindrical 
satellite is given by 

T = & x F (B 4) 

— Q, — 

where &_ is the distance vector from B* to the geometrical centroid 
(GC) of the satellite. 

In Appendix A, Figure A 1, the orbiting reference axes with unit 
vectors ' were defined The unit vectors n , 1 = 1,2,3, 

are defined m Section A of Chapter II (see Figure (2.l)). 

Since the GC is fixed m B, it is convenient to write 

A = 4^ + 4 2 n 2 -f 4^ (B.5) 
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where - £ n^, 1 = 1,2,3 If n is written m terms of n , 

i = 1,2,3, then with the use of (B 5) T will he m -terms of n , 

—a i 

l = 1,2,3, as desired. 

_ B*/a B* a B* . 

Since v = v - v = Vn and since v = rn-^ + r9nj, and 

a E • , . ^ ' 

v = x rn^ = ©^(ST x np^j where e £ is the angular speed of the 

earth and N is given in Appendix A, then 


^ = v^EL + 


r(0 - d E c s ) 




r °E S S C 0 


(B 6) 


where Sg, c^ are defined m Appendix A. 

From (B.6) it is seen that 

V 2 = f 2 + r 2 © 2 + r 2 [0 2 (c 2 + s 2 c 2 ) - 20^] 

With the aid of (C 7)-(C 9) equation (B 7) can he written as 

„S h 2 (l+e 2 + 2ecS) „ „ . d-d, 2 . 2 2, 

^ - - 2 7 7 2.2 L - 2h0 E°8 + r 0 E (c 6 + S 8 C 0 5 

a (1-e ) 

where c0 = cos 0 . 

If e ^ 0 05, then the expression for v differs from 


(B7) 


V = (h/a)(l + 2ec0) 1//2 = r0 


(B 8) 


only m the third significant figure. 

The substitution of (A. 5 ) into (B 6 ) and the results of this 
substitution into (B 2) produces an expression for the force, F, which 
is m terms n , 1 = 1,2,3. If this expression for the force is substi- 
tuted into (B with (B 5 ) and if the vector product is carried out, 

the result is T = T -.n + T ^n. + T „n„ where 

—a al— 1 a2— 2 a3— 3 


T _ = - 
al 


C D S P V 


{r(A 2 s 2 + ^ 3 c 2 s 3 ) - r(0 - 0 1J1 c R )[A o S-,c 


E 5 2T2 


+ ^ 3 (c 1 c 3 - S-^gSg)] + ^gCgSgti^Cg - ^(s-jCg + C;L s 2 s 3 )]} 

c D s P v 

2 *- r ^3 C 2 C 3 " £ 1 B 2^ + r ^ 0 " 0 E C S^' e 3^ S l S 2 C 3 + C l S 3^ +i l S l C 2^ 

, £- 0 ( Cont mued ) 


V 



+ r9 i‘W*s (B l’5 ' °1 S 2°3 ) - i l C l C 2 ]) 

C D SpV . 

a3 = 2 ^" r ^l c 2 s 3 ^2 C 2 C 3 ' + r ^ _0 E c &^f c l c 3 ~ S 1 S 2 S 3^ 

- e 2 U 1 s 2 o 3 + CjSj)] + re E o a s 5 [^ 1 (s 1 o 3 + o 1 s 2 s 3 ) 

- *2 ts l S 3 " =1 S 2 C 3 )15 (B 9) 


The quantities p and v m (B.9) are assumed to be suitably- 
given by (B.3) and (B 8) The quantities r,r, and 9 are given as 
functions of time m Appendix C The quantities C^, S, i = 1,2,3, 
pertain to a given satellite., and, the quantities h,a,Cg,s & ,0.p,r.p,K 
and pertain to a given orbit. Except for which is taken as 

= 2 2, given values of these quantities were chosen m Chapter II 
The angular speed of the earth is 9 = 7 3 x 10 ^ rad /sec. 

-fcj 

The error m the torque components, (B 9 ), are considered to be 

caused primarily by the error in p and the error m the product C^S. 

These errors were discussed above The total error in each of T , 

ai 

1 = 1 , 2 , 5 , is considered to be between 2000$ to large (night-time at 
500 miles) and 50$ too small (day-time at 5°0 miles) Generally., this 
means that the magnitudes of the components of T are between l/2 and 

“ 3 , 

20 times the correct values 


For the purpose of determining a feedback control law which will 

result m suitable attitude motion, when this motion is influenced by 

the atmosphere, this inaccuracy is insignificant Of course, performance 

measures such as fuel cost can be m error as much as T , 1 = 1,2,3, 

&.1 

are m error, but, generally, this depends on the magnitudes of the 
components of the other torques which influence the attitude motion. 
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APPENDIX C THE PARAMETERS OP THE ORBIT 
AS FUNCTIONS OF TIME 


King-Hele [23] and Sterne [33] have derived expressions for the 
rates of change of a , e, ~£!, 5 v and 0^ of a Kepler orbit due to the 
primary perturbing sources which are the earth's oblateness and the 
earth's atmosphere. For the orbits considered m this investigation 

-4 

the changes are less than 10 miles per day m a, 2 X 10 per day in 

e, 7° per day m £2, 1° per day m 6 and 10° per day m 9 For any 

3? 

one orbit these changes cause only an insignificant change m the terms 
m the attitude equations of motion Thus, since the equations of mo- 
tion are periodically time varying, they need to be integrated for a 

single orbit only so that for any one integration a, e, 0, 6 and $ 

P 

can be assumed constant. 

The angle 9 varies from zero to 2n radians For a Kepler or- 
bit 9 is related to r by 

r = (h 2 /p) (l + ec0) 1 (C l) 

- 2 
where p = m^G, c0 = cos0 and h is a constant such that r 9 = h 

Thus, 


^ = n 2 (1 + ee0) 3 (C 2) 

r J 


9 = n (l + ee0) 


2 


(C.3) 


and 



(l + ec6) 2 

2 / .. 2 \ 2 

a (1 - e ) 


(c.4) 


* For convenience c 9 = cos 9, but generally c0 = cos (0 + 0^). 
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£> o 2 2 

where n = p /hr and a (l - e ) = h /n 

The differentiation of C.3 "with respect to time and the substitu- 
tion of C 3 into the differentiated expression grves 

0 - -2n 2 es„ (l 4- ec0)"^ (C ♦ 5 ) 

u 

The angle 0 can be -written as an approximate function of t so 
that the right hand sides of equations C 1 - C 5 can be written as ap- 
proximate functions of t 

In the thesis e is assumed to be less than or equal to 0.03 In 
this case the function (l + ec0) can be expanded m a rapidly con- 
verging power series m e so that the integral of C 3 can be written 
as 



2ec0 + 


2 2 

3e c 0 - 


[REMAINDER]) d0 = n 



dt 


where (REMAINDER) = ij-e 3 [(c0 + ec 2 0) 3 /(l. + ee0) 8 ] _ for 0 < x < 0.05. 
If the maximum absolute -value of the REMAINDER is used, then 


0 


(REMAINDER) d0 i <6.5 


X 10 0 


0 


so what 


nt = 0(1-*- 3/2e 2 ) - 2e sm0 3/4-e 2 sm 20 
with an error less than 0 . 0 6 $ or 
nt = 0 - 2e sm0 


(C.6) 
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With an error only m the third significant figure (except at 9 = 0 
where the error is zero) 

Equations C 1 - C.5 are approximately given by 


r = (h 2 /|i) (1 - ec0) (C 7) 

^ = n 2 (1 + 3ec0) (C 8) 

r 

6 = n (1+ 2ec0) (C 9) 

r~ 2 = a -2 (l + 2ec0) (C 10) 

6 = -2n 2 es 6 (1 + 3ece) (C 11) 


with the approximation only in the third significant figure 
The use of C.6 and trigonometric identities results m 

cd = cos nt cos (2e sinfl) - sm nt sm (2e sine) (C 12) 


and 


s 9 = sm nt cos (2e sine) + cos nt sin (2e smfl) (C 13) 

The egressions for c0 and s9 given in C 12 and C.13 differ 
from exact values by less than OoOl. (This is easily seen by plotting 
[cos0 - cos (0 ± e)] and [sm0 - sm (0 ± e)] as a function of the 
error, e < 0 01 (e > 0), for various values of 6 ) Thus, if C 12 
and C 13 are substituted into C-7 - C 11, then C-7 - C 10 will be in 
error only m the fourth significant figure and C 11 will be m error 
by less than 1$ for most of the orbit. In C 7 - C 11 the factors 
(l - ec0), (l + 2ee0), (l + 3ec0) and s0 (l + 3ec0) can be written 
with the aid of C 12 and C 13 as 
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(C.l4) 


(l - ec0) = 1 - e cos nt cos (2e sinS) 

+ e sin nt sin (2e sm0) 

(l + 2ec0) = 1 + 2e cos nt cos (2e sm0) 

- 2e sxn nt sin (2e sin0) (C.15) 

(l + 3ec0) = 1 + 3e cos nt cos (2e sm0) 

- 3e sin nt sin (2e sm0) (C.l6) 

s0 (l + 3ec0) = sm nx cos (2e sxn0) + cos nt sm (2e sm0) 

2 2 

+ 3e cos nt sm nt [cos (2e sm0) - sm (2esm0)] 
2 2 

+ 3e (cos nt - sm nt) sm (2e sm0) cos (2e sm0) 

(C.17) 

Sj.nce e < 0.05, then |sm (2e sm0) | < 2e < 0.1 and |cos (2e 
sxn0) 1 = 1.00. Thus, C„l4 - C.l6 can differ from 

(l - ec0) = 1 - e cos nt 

(l + 2ec0) = 1 + 2e cos nt 

(l + 3ec0) = 1 + 3© cos nt 

only in the third significant figure and 

s0 (l + 3ec0) = sin nt + cos nt 0 2e 

2 2 

+ 3e (cos nt - sin nt) 

is m error by less than Vj>. Since sxn0 = sm nt + 2e sm0 cos nt 
(from C.13) "with less than lfo error for most of the orbit, then C.21 
can be -written as 


sin0 + 3 e cos nt sin nt 
• 2e sm0 (C.2l) 


(C.l8) 

(C.19) 

(C.20) 
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s0 (l + 3ec0) = sin nt + 5e cos nt sin nt 
or 

s0' < (l + 3ec0) = sm nt + 5/2 e sm 2nt (C.22.) 

with a total error of less than 1# for most of the orbit. (Actually, 
the difference m the exact values and the values given by C.22 is less 
than 0.01 for the entire orbit, but the percent error becomes unbounded 
as the exact value approaches zero), 

Now the substitution of C.l8 - C.20 and C.22 into C.7 - C.ll re- 
sults m 

r = r [l - e(eos nt -l)3 (C.23) 


- = n 2 (l + 3e cos nt) 
r^ 

0 = n (l + 2e cos nt) 

“2 -2 r \ 

r = a (1 t 2e cos nt) 


(C.24) 

(0.25) 

(C.26) 


o 

0 - -2n e(sm nt + 5/2e sm 2nt) 


(0.27) 


where m C.7 h /p has been replaced by a(l - e ) = r (l + e) and 
terms m e omitted. , The results m equations C.23 - C.26 differ 

from the exact values m the third or higher significant figure, and, 

• • 

the values for 0 which are given by C.27 are m error by less than 
OoOl parts in one for most of the orbit. 

If C.6 is used, then s0 and c0 can be written approximately as 
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(C.28) 


and 


sin (0 + 0 Q ) = sm (nt + 0 Q ) 

+ 2e cos (nt + 0 ) sxn nt (l + 2e cos nt) 

cos (0 + 0 ) = cos (nt + 0 ) 

- 2e sin (nt + 0 Q ) sin nt (l + 2e cos nt) (C.29) 

The order of approximation m C.28 and C.29 is the same as m C 12 

c.13. 
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appendix: d the state- space region of station-keeping 


The region of station-keeping, S, was defined m Section A of 
Chapter III Since it is clear that at least parts of the "boundary of 
S must he switching surfaces for a feedback control law, the choice of 
the region S from possible station-keeping regions should be a best 
choice. The criteria used for a best choice are small error, a minimum 
fuel expenditure and simplicity. 

If the angle, cp, between the local vertical and the satellite - 
fixed line [see Figure (2 l) ] and its rate, cp*, are required to 

be small, say M £ s, I cp' j sr, then the requirements on 0 , 1 = 2,3, 
and their rates are 





^ s 


W' 


* 


IVa + e s 8 3 1 


r 


+0 3 


lx 5 x 4 + X 5 X 6 1 < 


\[* 


r 


2 A 2 
X 3 x 5 


For near earth-point ing, the yaw angle and its rate must also be 
restricted to small values so that the requirements on 6 ^ and 0^ 
are | s s^, |@jj s s^. This region of station-keeping is not suf- 
ficient to avoid exceedingly large values m either 9 ^ or Also, 

the switching logic for a controller based on such a region is not 
simple. 

Generally, the station-keeping region should be closed region of 
the state-space which encloses the origin Thus, m any phase-plane 
projection of the region, the projection of its boundary should be a 
closed curve which encircles the origin Circles, ellipses, parts of 
parabolas, parts of straight lines, etc can be used for constructing 
the closed curves However, if parts of curves such as circles, para- 
bolas, etc., which are analytically described with squares of the state 
variables, are used, the error m the state variable signal from the 
sensors is compounded. Thus, the region which encloses the origin 
should be constructed with linear functions unless such a region results 


170 



m a significant fuel cost increase The region S, as given, is one 


such region, 
discussed m 


The dependence of the size of the region on fuel cost is 
Section B, Chapter VI 
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APPENDIX E. AN APPROXIMATE SOLUTION OF THE ADJOINT EQUATIONS 


For satellites which are "stable" in roll-yaw., the approximate 
solution of the adjoint equations for 



which is derived below is m excellent agreement with the digital 
computer solutions. 

Let t , t* ^ r s t*, denote the time at which the roll-yaw (pitch) 
e o e f 

trajectory re-enters the roll -yaw (pitch) projection of S if the 

trajectory has departed the projection of S Otherwise, let r = t*. 

e o 

With this fixed (after each encounter with the boundary of S) time, 
the piecewise linear equations, (4.4) and (4.6) with (4.7), can be 
approximated by piecewise constant equations by replacing t* with 
m A^, Ag and A 

Let t , ?* <. t <. ?*, i = 1, ...,6, Denote the time at which the 
1 o l f 

projected trajectory, x^(t*), exits the x^-projection of S. Then, 
since F(x i ) is lCP 1 , -10 ni or zero and since Al^, A^ and are 

now constants, equations (4.4) and (4.6) can be laplace transformed as 
follows 
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Pinch 


-F(x ) 

sP c - p c = — exp('-T ) - 3k A 

5 5e s 5s 336 

F(xg) 

^6 " p 6e = —T" e *P(- T 6 s) + p 5 


(E.2) 


where p = p (*r e ), 1 = 1,...,6, and P^ = P^s) is the Laplace 
transform of the variable p (t*)> 1 = 1 ,.,.,6, 

If equations (E l) are solved for P * 1 = 1> A* and equations 
(E.2) are solved for P , 1 = 5 , 6 , the results are 

P 1 - {p le s3 - (T 1 + p aPl + P 4e B 5 )s2 + [p le (B 4 ' V ' p 2e B 4 B 5 

- P 3e B 5 + VlfcW - B 1 T £ ■ B 5 T 4^ s + p le B 2 B 5 + P 2e^ B l B 3 " ^ r ) 

+ p 5e E lB2 f (B 3 - + B 5 T 3 + B 2 T t 

- + (B^B, - B |)T 2 - B 1 B 2 T 3 ]s‘ 1 )/A 

P 2 = t p 2e s3 + <P le - B 2 p 1k. + - < p 2e B 3 + p 3e B 2 + p 4e B 5 

+ T 1 + ® 2 V S - (p le B S + %e B 5 + B 3 T 2 ’ B 2 T 3 + B 5 V 
+ (B^ + B 5 T 3 )s" 1 }/A 

P 3 - lp 3e s3 + (B 5 P 2e " B 3 p lte ' + Cp le B 5 + P 2e B 3 B 4 

+ P 3e^ B l + B 2^ + P 4e B 5 B 2 + B 5 T 2 “ B 3\^ S + P le B 3 B 4 

+ p 3e B 4 B 5 + p l. e (E I B 3 + B 5> - Vl + Ws - (B 1 + B aV T 3 

- w. 1 ‘ [ Wi - Ws + < E f - iWV 6 ' 1 )/ 4 

P 4 - {p 4e s 3 + < p 2e B 4 + p 3e + + (p le B 4 + p 2e B 5 + p 4e B l 

+ \\ ' V s + p le B 5 + p 3eV' W * W * B A 

- (B 5 T 1 + B 1 T 3 )s“ ;L }/a (E 3) 
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P 5 ‘ (s 5e s - 


p 6e B 


’6 


“ V 


b 6 t 6 s 


■ 1 )/c 


s 2 + Eg) 


Pg - (P 6e s + P 5e + % - TjS-^/ts 2 + Bg) (S 4) 

]l O 

where - Ftx^expt-T^s), i = 1 ,...,6, A = s + (B^ - ® 3 + BgB^)s 

+ B 5 CB 2 - V s + ®5 - B 1 B 3 “ d B 1 - k l4 B 2 - *2*1’ B 3 ■ k 2<4 + 3A 5 >’ 

= K^A^; B^ = -Ag and Bg = Sk^Ag with A^, i = 1 , 2 , 5 , evaluated 
at t* = t . 

e 2 
Let B^ - = 2a, Bg - B^ = b and B^ - B^B^ = c. Then 

A can he written as 


A = s + 2as + 


B.-bs 

5 


+ e 


(E.5) 


If the characteristic equation, (E. 5), has roots which do not change 

type, e.g. ; from type imaginary to type complex, with time, then equations 

(E.3) and (E.4) can be inverted. once -and-for-all to give the approximate 

solution of (4.4) and (4.6) for the entire time of one orbit. This is 

the case if B,_b is zero or very small. 

For satellites (l) and (2) B_b varies periodically from about ELCT^ to 
_1|_ y 

-10 and the roots of (E.5) have real parts less than 0.1 m magnitude 
Thus, the exponential factors in the solution due to these real parts 
are nearly unity for solution times i.e., the time intervals between two 
successive encounters with the boundary of S, less than unity (in 
Chapter IV it is seen that dozens of encounters occur'- m the time of 
one orbit, which takes about six time units to complete.) 

If the term m (E.5) with B^b as a factor is omitted, then (E.5) 

2 2 2 2 i- 2 1 1 

can be wri tten as A = (s + oi)(s + (xt) where ax = a + v a ~c and 

2 / ^ ^ ^ 
o> 2 = a - va -c . If this expression of A is substituted into equations 

(E.3) and if, on inverting (E„3) and (E.4) for positive and real values 

of and Bg, the functions sm y and cos y are replaced by 

y and 1 - y^, respectively, then the approximate solution of equations 

(4.4) and (4.6) is 
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P x = p le Cl -[(B 3 - B 4 ) + (o> 2 + ai 2 )]^*) 2 } _ P 2e [B 1 AT* 

+ B^At*) 2 ] - P 3e B 5 (AT^) 2 - p 4e [B 5 AT* - B^At*) 2 ] 

- F(xj^)B 5 (at ||;) 2 - F(x 1 )At£ - FU^B-jAt * 2 

P 2 = P le AT * + ^e 1 - 1 + ^ + b 3 Hat*) 2 } - P 3e B p (AT*) 2 

- P 4e ( B 2 A^ + B^t* 2 ) + F(x 2 )At* - F( X;l )(At *) 2 - F(x 4 )B 2 (At |^) 2 

P 3 = p ie B 5^ AT ^^ 2 + P 2e*- B 5 ATJf + B 3 B 4 (At*) 2 ] + ^Se^ 1 ^ 

+ (B 1 + B 2 )](At*) 2 } - p^ e (B 3 AT^ + B 2 B 5 (At *) 2 + F(x 2 )B 5 (At ^) 2 

- F(x 3 )Atj - F(x^)B 3 (At ^) 2 


P4 = p ie B 4^ AT ^^ + p 2e [B 4 ATJf + + p 3e ATJf 

+ p 4 e £ 1 "[(“l + “2 ^ + V^*) 2 } + F(x^)At|^ - F(x 3 )(At^) 2 
+ F(x 2 )B^(AT£) 2 (E 6) 

P 5 = V 5e^ 1 “ ~ P 6e B 6 (AT*) 2 “ f ( x 5 )A t J - F(x 6 )Bg(AT^) 2 

P 6 = P 5e AT * + Fge'- 1 “ B 6( At *) 2 -1 + F ( X 6^ AT ^ “ F^K^t*) 2 (E 7) 

where At* = t* - t^ and At* = t* - t^, 1 = 1, ,6 

An approximate solution of the adjoint equations for "unstable" 
satellites can be derived m a similar manner However, instead of sine 
and cosine functions (which can be replaced by simplier functions) 
appearing m the equations, the equations contain the exponential function 
and are much more complicated m yaw-roll In Section A, Part 3 of 
Chapter IV an approximate solution is presented for "unstable" pitch motion 
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APPENDIX F. DIGITAL AND ANALOG COMPUTER PROGRAMS 


In Figures P 1, F 2 } F 5 and F 6 are given listings of the digital 
computer programs. These are given m the order m -which they are used 
m the text The word "clock" which appears m the digital computer 
program, is an ALGOL procedure This procedure was used to determine the 
elapsed time required for the execution of certain parts of the programs 
so that measures could he taken to reduce excessive computation times 
The symbols used m the differential equations (DE) procedures are not 
the same as m the differential equations m the text since these symbols 
were reserved for the plot routine. However, the equivalences are given 
m the comments 

In Figures F 3 and F 4 are given simplified analog computer programs 
which were used to simulate the acquisition motions from S + to S and 
the steady-state motion for the suboptimal steady-state control. The 
time delays required m these simulations were obtained from the time 
delays m the comparators by time scaling 
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PROCEDURE KUTTAMERSON(N» X#HH# Y # F# EPS# AB # ERROR# STEPSIZE )! 

VALUE N#HH# EPS#AB # STEP SI ZEHNTEGER Ni REAL X#NH# EPS#AB ! 

REAL ARRAY YC03! PROCEDURE Ft BOOLEAN ERROR# STEPSIZE! 

COMMENT! VERSION OF 660518 660722 

EPS AND AB ARE THE RELATIVE AND ABSOLUTE ERROR BOUNDS RESP. 
STEPSIZE TRUE TO WRITE STEPSIZE WHEN CHANGED 
SIEPSIZE FALSE FOR NO OUTPUT 

ERROR IS SET TRUE IF STEPSIZE BECOMES TOO SMALL ELSE FALSE! 
BEGIN COMMENT KUTTA MERSON INTEGRATES A SYSTEM OF N FIRST ORDER 

ORDINARY DIFFERENTIAL EQUATIONS. SEE L. FOX# "NUMERICAL 
SOLUTION OF ORDINARY AND PARTIAL DIFFERENTIAL 
EQUATIONS"# P, 24# PERMAGON PRESS# 1962 ! ' 

OWN REAL HC# FINAL# H2# H3, H6# H8# ERR# TEST# T, Hi 

OWN INTEGER I#CU#CUTi OWN BOOLEAN DBLI LABEL L# KM# RETURN! 

OWN REAL ARRAY Yl# Y2# FO# FI# F2E0*303! 

COMMENT EXCEPT FOR HC# THE OWN VARIABLES ARE FOR SPEED ONLY! 

FORMAT MSSG("TH£ STEP SIZE IS NOW"# R12.5#" AT T="#R12.5)i 
DEFINE FORI = FOR 1*1 STEP 1 UNTIL N DO f# 

CONSTANTS = H2+H/2.0! H3«-H/3.0! H6*H/6,0i H8«-H/8.0 #i 
COMMENT CHECK FOR INITIAL ENTRY AND ADJUST H IF NECESSARY ! 

ERROR * FALSE! 

H «■ HH ! 

IF N=0 THEN BEGIN HC * Hi GO TO RETURN END! 

IF H=0 THEN GO TO RETURN! FINAL «• X+Hi 
IF HC=0 THEN HC * H! 

IF EPS*0 AND ABS(H)>ABS(HC) THEN 

IF SIGNCH)^SIGN(HC) THEN H <- HC * -HC ELSE H * HC! 

COMMENT! CUT IS THE NUMBER OF TIMES THAT THE STEP SIZE IS ALLOWED TO 
HALVE ITSELF IN SUCCESION! 

CUT * 4! 

CU «• CUT! 

T «• X + H! X * FINAL! CONSTANTS! 

COMMENT MAIN KUTTA-MERSON STEP LOOP ! 

L!FGR T«-I STEP H UNTIL FINAL DO 
BEGIN KM! F(T-H#Y#FO>i 

FORI Yim * FOCI3XH3+YEI3! F(T-2xH3# Yl# Fl)i 

FORI Yim * (F0m + FltI])xH6 + YtI]i F(T-2*H3, Yl# Fl)i 

FORI Yim * <Fin]x3,0+F0m)xH8 + Ymi FCT-H2# Yl# F2)i 

FORI Yim * CF2mx4.0-Fimx3.0+FOm>XH2+Ytm F(T# Yl# FI)! 

FORI Y2CI] * (F2mx4,0+Fim + F0[I))XH6+YEIJ! 

COMMENT DUES THE STEP SIZE H NEED TO BE CHANGED ! 

IF tPS*0 THEN 
BEGIN DBL * TRUE! 

FORI BEGIN ERR*ABS(Yim-Y2m)xO,2! TESTfABSCYl E 13 JxEPS! 

IF ERR>TEST AND ERR>AB THEN COMMENT HAlF H! 

BEGIN H * H2! T *T -H2 ! 

IF CCU*CU-1)<0 THEN BEGIN ERROR * TRUE! EPS*0! 

GO TO KM! END! 

IF STEPSIZE THEN WRITE ( MSSG# H# T ) ! 

IF T+H = T THEN BEGIN X*Ti ERROR «• TRUE! GO TO RETURN! 
ENO! 

CONSTANTS! GO TO KM! 

ENO# 

IF 64 * OxERR>TEST THEN DBL * FALSE! 

END! 


Figure F 1. Solution of Approximate# Minimum-Fuel Optimal# Station- 
Keeping Equations [(4 3) and (4 4 )] for Roll -Yaw. 
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PROCEDURE DECTAU.Y.DY)! 

VALUE TAUI 
HEAL TAU! 

ARRAY Y,DYI03! 

BEGIN 

REAL Ai*A2>A3»A4! 

COMMENT YC13=1000xXl* Y[23=lOOOxX2, Y[ 3]= 1000 X X3» Y [ 4 3 = 1 OOOXX4 > 
Y[53 = O.OixPi> YC63 = 0 ,OlxP2/ Y [ 7 3 = 0 , 0 IxP 3 > Y 1 8 3 =*0 . 0 lxP H, 
Y[93=J,R0LL-YAWsjRY! 

Yl* YtUJ Y2«- Y [ 23 1 Y3«- YE33! Y4 + Y[43! 

Y5+- YtbJ! Y6 + YI63; Y7 + YC73! Y8«- Y[83! 

RG<-Q , If 

FXH IF A8S(Y1)SRG THEN 0.0 ELSE IF Y1>RG THEN PR1 

ELSE "PR 1 1 

FX2«- IF ABS( Y23SRG THEN 0.0 ELSE IF Y2>RG THEN PR2 

ELSE -PR2J 

FX 34- IF ABSCY33SRG THEN 0.0 ELSE IF Y3>RG THEN PR3 

ELSE -PR3J 

FX4«- IF ABS(Y4)<RG THEN 0.0 ELSE IF Y4>RG THEN PR4 

ELSE -RR4! 

COMMENT W 1 = V 1 > W2=V2! 

Wl* IF ABS(Y6)>0»01 THEN NixSIGN(Y6) ELSE 0,01 
W2«- IF A8S(Y8)>0.01 THEN N2XSIGNC Y8 ) ELSE 0.0! 

CT«-COS( TAU)! ST + SINCTAU )! 

A1«-1 + 4xExCT, A24.2xEx(ST+5x£xSTxcT); 

A3*1 + 2xExCTI A4*4+13xExCT! 

DY [ 1 ]<• -Y2! 

DYC2 34- KlxAlxYl-A2xY3-K3xA3xY4-1000xwiI 
DYI33* -Y4! 

D Y E 4 3 *• A2xYi + K4xA3xY2-K2xA4xY3"1000xW2j 
DYC53+ -0,01xFXl-KlxAixY6-A2xY8; 

DY 1 6 3*- 0,0IxFX2 + Y5-K4xA3xY8i 
0 Y 1 7 H -0,0lxFX3+A2xY6 + K2XA4xY8I 
DYC83* 0,01xFX4+K3xA3xY6+Y7; 

DYE93* ABSCHl )+ABS(W2)! 

END DE> 

COMMENT INITIAL(FINAL) CONDITIONS! 

START » READ{E # Ni#N2>Kl,K2,K3,K4,PRl.PR2^PR3/PR4.X10>X20.X30,X40. 

MUHFINISH3! 

TAU«-TEOJ«-0! 

Ym + ioooxxioj xito3<-xio; 

Y E 23 «■ 1000xX2Q! X2E03*X20! 

YE334-1000XX30! X3[0]*X30! 

Y [ 4 3 «■ 1 OOOXX4Q! X4[03*X40! 

YC53* IF ABSCX10)=RG THEN -MUxSIGN< XlO) ELSE 0.0! 

Y E 6 3 «• IF ABS(X20)=RG THEN -MljxS I GN ( X20 > ELSE 0.0! 

P2C03<-Yt6]! Pi[0l«-Yt53! 

Y 1 7 3 *• IF ABSCX30)=RG THEN -MUxS I GN ( X 30 ) ELSE 0.0! 

YC8J* IF ABS(X40)=RG THEN -MUxS I GN C X40 ) ELSE O.j)! 

P4t03<-YC83! P3E03«-YE71! 

Y 1 9 3 *-0! 

WRITECLA8L), 

HRITE(RESL>TAU>X10/X20>X30>X40>Yr63>YE8],Wl>H2)! 

COMMENT CALCULATING USING KUTTAMERSON AND LOADING ARRAYS! 

FOR L* 1 STEP I UNTIL 628 DQ BEGIN 


Figure F . 1 Continued. 
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KUTTAMERS0N<9#TAU#0.0l#Y>0E>e-4#e«*5#FINISH#FALSE); 

X1CL]«-0*001xYCI 31 X2[L3*0,001xYC 23; X3 t L 3«-0 . 0 0 1 x y [ 3 3 } 
X4[U<-0.001xY[43I 

P2CLI«-100xY[ 63; p4[L]«-100xYt8]J P 1 CL ]<■ lOOxYI 5 ] f P3CL 3*100xYC7 JJ 
PRINT<-1I 

IF L MUD PRINT = 0 THEN 

WRlTEtRESl^T[L]#Ytn#Y[2]#YC33#Y[43#YE63#YC8]#Wl,W2)I 

END KUTTAMERSON LOOP! 

JR Y*Y I 9 j I 

WRITE(PARM#E#K1#K2#K3#K4> JRY); 
WRITECtPAGE3#PARMT#PRl#PR2#PR3#PR4#Nl#N2#MU>; 


Figure F 1 Continued. 
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PROCEDURE KUTTAMERSONCN* X*HH* Y* F, EPS* AB * ERRORS STEPSIZE); 

VALUE N*HH* EPS,AB * STEPS I Ze; INTE GER N1 REAL X*HH* EPS*AB ) 

REAL ARRAY YC03; PROCEDURE FJ BOOLEAN ERROR* STEPSIZE; 

COMMENT! VERSION OF 660518 660722 

EPS AND A8 ARE THE RELATIVE AND ABSOLUTE ERROR BOUNDS RESP. 
SIEPSIZE TRUE TO WRITE STEPSIZE WHEN CHANGED 
STLPSIZE FALSE FOR NO OUTPUT 

ERROR IS SET TRUE IF STEPSIZE BECOMES TOO SMALL ELSE FALSE! 
BEGIN COMMENT KUTTA MERSON INTEGRATES A SYSTEM OF N FIRST ORDER 

ORDINARY DIFFERENTIAL EQUATIONS, SEE L* FOX* "NUMERICAL 
SOLUTION OF ORDINARY AND PARTIAL DIFFERENTIAL 
EQUATIONS"* P, 24, PERMAGON PRESS* 1962 t 
OWN REAL HC, FINAL* H2* H3* H6* H8* ERR* TEST* T* H* 

OWN INtEGLR I*CU*CUTJ OWN BOOLEAN DBLJ LABEL L* KM* RETURNS 
OWN REAL ARRAY Yl* Y2* FO* FI* F2COt30]J 
COMMENT EXCEPT FOR HC* THE OWN VARIABLES ARE FOR SPEED ONLY } 

FORMAT MSSG ( "THE STEP SIZE IS NOW"* R12.5*" AT T="*RI2,5)F 
DEFINE FORI = FOR I<-1 STEP 1 UNTIL N DO rf* 

CONSTANTS = H2«-H/2,0; H3«-H/3.0; H6«-H/6.0* H8«-H/8.0 *} 

COMMENT CHECK FOR INITIAL ENTRY AND ADJUST H IF NECESSARY * 

ERROR «■ FALSE! 

H «- HH / 

IF N=0 THEN BEGIN HC * H* GO TO RETURN END; 

IF H = 0 THEN GO TO RETURN I FINAL <• X + HT 

IF HC = 0 THEN HC «• HI 

IF EPS?0 AND ABSCH)>ABSCHC) THEN 

IF SIGNCHXSIGNCHC) THEN H «■ HC * -HC ELSE H <■ HC; 

COMMENT! CUT IS THE NUMBER OF TIMES THAT THE STEP SIZE IS ALLOWED TO 
HALVE ITSELF IN SUCCESION! 

CUT <• io; 
cu *• cut; 

t *■ x+h; x «• final; constants; 
comment main kutta-merson step loop ; 

L ! FOR T*T STEP H UNTIL FINAL DO 
BEGIN KM I F(T-H,Y*FO); 

FOR! YHI3 <• F0CI]xH 3 + YCI]; FCT-2xH3* Yl* FD; 

FORI Y1CIJ «■ (FOCIJ + Fim)xH6 + YCI3; FIT-2XH3* Yl* FI); 

FORI Y1CI3 «• (Fimx3,0 + F0m )XH8 + Yin; FCT-H2* Yl, F2); 

FORI Yim <• CF2CI3X4,0-Fimx3,0+F0[I])XH2 + Ytn; FCT* Yl* Fl>; 
FORI Y2CIJ «■ (F2mxii,0 + Fim + F0[I3)XH6 + Ytn; 

COMMENT DUES THE STEP SIZE H NEED TO BE CHANGED ; 

IF LPS*0 THEN 
BEGIN OBL <• TRUE; 

FORI BEGIN ERR<-ABS(Yim-Y2II])x0.2; TEST* ABS C Yl [ 1 3 ) xEPSI 
IF ERR>TEST AND ERR>AB THEN COMMENT HALF H; 

BEGIN H 4- H2; T«-T-H2; 

IF CCU4-CU-1X0 THEN ERROR «■ TRUE; 

if stepsize then writecmssg*h*t); 

if t+h=t then begin x*-t; error «• true; go to return; 

end; 

constants; go to km; 
end; 

IF 64 • OxERR>TEST THEN DBL ♦- FALSE; 
lnd; 

IF DBL AND H < HH THEN BEGIN H <- 2,0xHj 


Figure F 2 Solution of Approximate, Minimum-Fuel Optimal, Station- 
Keeping Equations [(4*5) and (46)] for Pitch 
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IF STEPSIZE THEN WRIT£CM5SG»H>T); 

cu <■ cut; 

CONSTANTS END DOUBLE HJ 

end; 

FORI YCI] «• Y2tli; 

end kutta merson loop; 

IF £PS=0 THEN GO TO RETURN; 

comment now be sure to have t = final ; 
hc <■ h; h * final-(t-H); 

IF ABSCH)>ABS(FINAL)xi,A55l9l5228gi-li THEN 

begin t «• final; eps <■ o; constants; go to l end; 

RETURN! ENU KUTTA MERSON; 

PROCEDURE D£CTAU*Y»DY); 

value tau; real tau; 

ARRAY Y*DY{0]; 

BEGIN 
REAL A*8; 

CT*COSCTAU); ST«-SIN(TAU>; 

A<-2xEx(ST + 5xExSTxCT); 
b«-i+3xe*ct; 

COMMENT YU] = 1000XX5^ YI2] = 1000xx6> Y l 3 ] «0 . 0 lxP5, Y [ 4 ] *0 .0 1 xP6 , Y C 5 3 = 
j,pitch=jp; 

Y i «■ YC13, Y2* Y[2j; Y3* Y[3]T Y4* Y[43; 

RG5<-0.U RG6*0«i; 

FX5+- IF ABS( Y 1 )<RG5 THEN 0 ELSE 

IF Y 1 >RG5 THEN PR5 ELSE -PR5; 

FX6* IF AB$(Y2)<RG6 THEN 0 ELSE 

IF Y2>RG6 THEN PR6 ELSE -PR6I 

COMMENT W3=V3; 

W3+ IF ABSCY4)>O.Oi THEN N3xSIGN(Y4) ELSE 0; 

DYI 1 ]«• -Y2I 

D Y [ 2 3 <• 3xK3xBxyi + ioOOX{A"W3); 

DY C 3 ]<■ -0.0lxFX5-3xK3xBxY4; 

DYI4]«-0.01xfX6 + Y3; 

DY[5 ]«■ ABS( W3 ) ; 

END DE > 

comment initialcfinau conditions; 

START! READ<E,K3,N3,PR5,PR6,X50>X60* YC3I,YC4]>YC53 > [FINISH]; 

tau*o,o; uo3*o.o; 

YI nuoooxxso; 

X5I03«-X50; 

Y[2]«-1000XX60; 

X6I0]«-X60; 

YC33<- IF ABSC X5Q } = RG5 THEN “MUxSIGN(X50 ) ELSE 0.0; 

Y C 4 3 *■ IF ABSC X60 ) = RG6 THEN -MUXS I GN C X60 ) ELSE 0.01 
P5[0]*YC33; 

P6[Q]*YU3; 

WRITECLABL); 

WRITECHESL* TAU*X50,X60, Yt33>Yt43* W3); 

COMMENT CALCULATING WITH KUTTAMERSON AND LOADING ARRAYS; 

FOR L * 1 STEP 1 UNTIL 620 DO BEGIN 

KUTTAMERSQNC5UAU,0,01*Y>OE>e-4>e-5> NUTS# FALSE); 

IF NUTS THEN 

WRITEC<"STEPSIZE WAS CUT FOUR TIMES AND KMC CGNTINUED>T="> 
F6,3*>>TA0); 


Figure F.2. Continued. 


l8l 



X5CL-H-0.001XYU3J X6 t LI *0 , 00 lxY t 2] I 
P5[LJ«-100*YE3]> P6IL]«-100xY[4]I 
T[L3«-Li 
PRINT+U 

IF L MOD PRINT = 0 THEN 
HRnE(RESL,TCL3,Yin*Yt2I>YC3],Ym>W3); 

END CALCULATING and LOAOING LOOP! 

JP«-H5 U 

WRITEUP AGE], PARM,E*K3>N3,PR5>PR6>MU.»JP); 


Figure F 2 Continued. 
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Figure F A Analog Computer Program for Pitch (Generator of Aerodynamic 
Torque Term Hot Shown) 
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PROCEDURE KUTTAMERSONCN, X, H , Y, E, EPS, AB , ERROR, STEPSIZE); 

VALUE H, EPS, A8 , STEP SI2EJ INTEGER N } REAL X, H, EPS,AB ; 

REAL ARRAY YlOi; PROCEDURE F; LABEL ERROR; BOOLEAN STEPSIZEJ 
COMMENT EPS AND AB ARE THE RELATIVE AND ABSOLUTE ERROR BOUNDS RESP. 
STEPSIZE TRUE To WRITE STEPSIZE WHEN CHANGED, FALSE FOR NO OUTPUT I 
BEGIN COMMENT KUTTA MERSON INTEGRATES A SYSTEM OF N FIRST ORDER 

ORDINARY DIFFERENTIAL EQUATIONS. SEE L. FOX, "NUMERICAL 
SOLUTION OF ORDINARY AND PARTIAL DIFFERENTIAL 
EQUATIONS", P, 24, PERMAGON PRESS, 1962 f 
OWN REAL HC, FINAL, H2, H3, H6, HB, ERR, TEST, T; 

OWN INTEGER IJ OWN BOOLEAN D8L> LABEL L, KM, RETURN, 

OWN REAL ARRAY Yl, Y2» EO, El, E2C0J30]; 

COMMENT EXCEPT FOR HC, THE OWN VARIABLES ARE FOR SPEED ONLY! 

FORMAT MSSG("THE STEP SIZE IS NOW", E18.ll); 

DEFINE FORI = FOR 1*1 STEP 1 UNTIL N DO #, 

CONSTANTS = H2*H/2.0; H3*H/3,0; H6*H/6.0; H8*H/0.O #; 

COMMENT CHECK FOR INITIAL ENTRY AND ADJUST H IF NECESSARY ; 

IF N=0 THEN BEGIN HC * HI GO TO RETURN END; 

IF H=0 THEN GO TO RETURN ; FINAL * X+HS 
IF HC=D THEN HC * h; 

IF EPS*Q AND ABS(H)>A8S(HC) THEN 

IF SIGN(H)^SIGNCHC) THEN H * HC * -HC ELSE H * HOT 

t * x+h, x * final; constants; 

COMMENT MAIN KUTTA-MERSON STEP LOOP ; 

L J FOR T*T STEP H UNTIL FINAL DO 
BEGIN KM J ECT“H,Y,FO); 

FORI YltlT * FOlI]xH3 + Ym; FCT-2xH3, Yl, FI); 

FORI YltlT * CF0U]+FlCI])xH6+Ym; FCT-2XH3, Yl, FI), 

FORI Y1CI1 ♦ (Fimx3.0 + F0m)xH8 + Ym; FCT-H2, Yl, F2); 

FORI YUIT * (F2mxzi,0-FUnx3.0 + Fon])xH2 + YCl); FCT, Yl, FI)) 
FORI Y2 1 1 3 * CF2mx4,0*Fim + FOm)xH6+Y[I3, 

COMMENT DOES THE STEP SIZE H NEED TO BE CHANGED ) 

IF EPS^O THEN 
BEGIN DBL * TRUE) 

FORI BEGIN ERR*ABS(Yim-Y2tIJ)x0.2; TEST* ABS C Y 1 [13 )x£P SI 
IF ERR>TEST AND ERR> AB THEN COMMENT HALE H) 

BEGIN H * H2) T*T-H2; 

IE STEPSIZE THEN WRI TEC t DBL3 , MSSG, H)3 
IF T+HsT THEN BEGIN X*TI GO TO ERROR £N0> 

constants; go to km; 
end; 

IF 64,0xERR>TEST THEN DBL * FALSE; 

end; 

IF DBL THEN BEGIN H * 2xh; 

IF STEPSIZE THEN WRITE C C DBL 3 , MSSG, H); 

CONSTANTS END DOUBLE H; 

end; 

EQKI YII3 * Y2II3; 

END KUTTA MERSDN LOOP; 

IE EPS = 0 THEN GO TO RETURN; 

COMMENT NOW BE SURE TO HAVE T = FINAL ; 

HC * HI H * FINAL-CT-H); 

IF ABSCH)>ABSCFINAL)xt.455l915228P-H THEN 

BEGIN T * FINAL! EPS * 0, CONSTANTS; GQ TO L END; 

RETURNi END KUTTA MERSON; 


Solution of Nonlinear, Suboptimal Acquisition Equations 
[(2 1*0~(2 16 ) and (5 2)] 
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IF DBL and H < HH THEN BEGIN H * 2.0XHJ 
IF STEPSIZE THEN WRITECMSSG/H/T)! 

CU * CUT ! 

CONSTANTS END DOUBLE HI 
END* 

FORI YCI] * Y2CI3! 

end kutta mersun loop! 

IF EPS S Q THEN GO TO RETURN/ 

COMMENT NOW BE SURE TO HAVE T = FINAL ! 

HC * H! H * FlNAL-CT-H)! 

IF ABS(H)>ABSCFINAL)xl,A55l915228B-ll THEN 

BEGIN T * FINAL! EPS * 0! CONSTANTS! GO TO L END! 

RETURN! END KUTTA MERSQN! 

PROCEDURE DE(TAU/Y/DY}! 

VALUE TAU! REAL TAU! 

ARRAY Y> D Y 1 0 1! 

begin 

REAL A1/A2/A3/A4/TA/C1/C2/C3/S1/S2/S3/CT/ST! 

T A* TAU + THETAO! CT* COSCTA)! ST* SINCTAJI 

A 1 * l + 2x£xCT! A2*-2x£xST! A3*l+3xExCT! A4* CxAlxEXP(KxExCCT-l ) )! 

COMMENT YUJsXl/ Y[23=X2/ Y[33=X3/ Y[4l=X4# Yt5]=X5, Y C 6 ]=X6# Y[71=J! 
Yl*Ytl3!Y<i*Yt23! Y3 + Y [ 3 3 ! Y4*Y[43J Y5*Y[53! Y6*Yt63! 

C1*C0SCY13! C2*C0S(Y3)! C3*COSCY53! 

SH-SINCY13! S2*SIN<Y33! S3*SINCY53! 

VI* IF ABSI Y 1 3 <0 *0002 AND ABSC Y2X0.0002 THEN 0 ELSE 

• CNl/2)x(SlGNCY2XABSCY23 + 20xYl) + SIGNCY2 + MlxYlxSlGN(Qi-ABSCYl ))) 

V2* IF ABSI Y3)<0, 0002 AND ABS C Y4 ) <0 . 0002 THEN 0 ELSE 

-{N2/2)x(SIGN(Y4xaBSCYA ) + 20xY3) + SIGNC Y4 + M2xY3xSlGN(Q2-ABSCY3) )3 >! 
V3* IF ABSCYi>)<0,0002 AND ABSC Y6 )<0. 0002 THEN 0 ELSE 

-(N3/2)x(SlGN(Y6XABSCY6)+20xY53+SIGNCY6+M3xY5xSlGNCQ3"ABSCY53))3! 
Wl* Y2xC2xC3+Y4xS3+Alx(SlxS3-ClxS2xC33! 

H2* Y4xC3-Y2xC2xS3+Alx(SlxC3+ClxS2xS33! 

W3* Y6+Y2XS2+A1XC1XC2! 

DYE13+ Y2! 

DYE23*(Y2xY4xS2-Y4xY6+Alx( Y4xClxC2-Y2xSlxS2-Y6xsi )+A2xcixS2+K2xS3xW3xWl 
-K1 xC3xW2xW3-3xa3x(K1+K2)xC2xs2xC3xS3+A4xc3x( S 1XC2+C1XC3) 
+C3 x Vl"S3xV2)/C2! 

OY [ 3 3* Y4! 

D Y [ 4 3 * Y2XY6xC2-AlxClx(Y2+Y6xS2)-A2xsi-w3x{KlxW2xS3+K2xWlxc3) 
+3xA3xCK2xC3xC3-KlxS3xS3)xC2xS2+C3xV2-S3xVi! 

DYI53* Y6! 

DY[63*-Y2xT4xC2-DYC23xS2+AlxCY2xSlxc2+Y4xcixS2)-K3XHixW3"3xK3xA3xC2xC2x 

C3xS3“A2xClxC2+A4xClxCS3-C3)+V3! 

DY[73* ABS(V1)+ABSCV2)+ABS(V33! 

END OE! 

COMMENT INITIAL CONDITIONS! CLOCK! 

START! REAUCE/Kl # K2/K3/Nl/N2/N3/X10/X20»X30>X40/X50/X60/Ml/Qi/M2/Q2/ 

M3/S3/C/0>K,THETA0/XlL/X2L/X3L/X4L/X5L/X6L/XlRj.X2R/X3R/X4R/X5R»X6R) 

[FINISH]! 


TAU* Y 1 

[7 3* II 

'03* 0] 

X 1 C 0 3 * 

YC1J* 

X10! 

X2C03* 

Y[23* 

X20! 

X3 [ 03* 

Y[33* 

X30! 

X4C03* 

Y[4J* 

X40! 

X5 [ 0 3* 

YC53* 

X50! 


Figure F 5 Continued. 
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X 6 C 0 3 Y[6J«- X60! 

WRITECLABL3! 

WRITE(RESL*TAU*X10*X20.»X30.»X40#X50*X60)J 

COMMENT CALCULATING WITH KUTTAMERSGN AND LOADING ARRAYS! 

FOR L* 1 STEP 1 WHILE LSIOOO AND CABSCYt23)>D OR ABSCY[4])>D OR 

ABSC YI6 ] 3>D) DO 

BEGIN 

KUTTAMERSUN(7*TAU>0.002> Y»DE*B- Up 0-5* NUTS* FALSE)! 

IF NUTS THEN WRITE(<"STEPSIZE WAS CUT FOUR TIMES BUT KMC CONTINUED* 
T aW *F6.3*>*TAU)! 

X1CL3* Y[13! 

X2CL3* Y 1 2 1 ! 

X3EL3* Y E 3 I ! 

XAEL3+ Y C A 3 ! 

X5CL3* Y I S 3 ! 

X6EL3* YC63! 

TCL3«-0,002xl! 

PRINT <-1 ! 

IF L MOD PRINT=0 THEN HR ITEC RE SL» T CL 3 * Y [ 1 3 * Y [ 2 3 * Y t 33 * Y [ 4 3> Y I 5 3 * Y 1 6 3 3 ! 

END CALCULATIND AND LOADING LOOP! 

U* Y C 7 3 ! 

WRITE! [PAGE3*PARM*E*K1*K2»K3*C*D*J3! CLOCK! 


Figure P.5. Continued 
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PROCEDURE MJTTAMERSQNCN> X#HH# Y * F # EPS# AB * ERROR# STEPSIZE3# 

VALUE N*HH> IPS#A<3 * STEPSIZE* INTEGER NJ REAL X*HH* EpS#AB * 

REAL ARRAY Y C 0 I # PROCEDURE F # BOOLEAN ERROR# STEPSIZE# 

COMMENT. VERSION OF 660518 660722 

EPS AND AB ARE THE RELATIVE AND ABSOLUTE ERROR BOUNDS RESP, 
STEPSIZE TRUE TO WRITE STEPSIZE WHEN CHANGED 
SltPSIZE FALSE FOR NO OUTPUT 

ERROR IS SET TRUE IF STEPSIZE BECOMES TOO SMALL ELSE FALSE! 
BEGIN COMMENT KUTTA MERSON INTEGRATES A SYSTEM OF N FIRST ORDER 

ORDINARY DIFFERENTIAL EQUATIONS. SEE L. FOX# "NUMERICAL 
SOLUTION OF ORDINARY AND PARTIAL DIFFERENTIAL 
EQUATIONS"# P. 24# PERMAGON PRESS# 1962 # 

OWN REAL HC# FINAL# H2# H3# H6# H8# ERR# TEST* T# H* 

OWN INTEGER I#CU#CUT# OWN BOOLEAN OBL# LABEL L» KM# RETURN# 

OWN REAL ARRAY Yl# Y?# FO# FI# F2C0.303! 

COMMENT EXCEPT FOR HC# THE OWN VARIABLES ARE FOR SPEED ONLY! 

FORMAT MSSG{"THE STEP SIZE IS NOW"# R12.5#" AT T="*R12.5)! 

DEFINE FOHI = FOR 1*1 STEP 1 UNTIL N DO ** 

CONSTANTS = H2*H/2.0# H3*H/3.0# H6*H/6.0# H8*H/8,0 t, 

COMMENT CHECK FOR INITIAL ENTRY AND ADJUST H IF NECESSARY 1 
ERROR * FAl SE# 

H «- HH # 

IF N=0 THEN BEGIN HC * H# GO TO RETURN END# 

IF H=0 THEN GO TO RETURN! FINAL <• X + H# 

IF HC=U THEN HC * H# 

IF EPS^O AND AB S { H 3 > ABS ( HC ) THEN 

IF SIGNCH)#5IGN(HC) THEN H <- HC * -HC ELSE H «- HC! 

COMMENT! CUT IS THE NUMBER OF TIMES THAT THE STEP SIZE IS ALLOWED TO 
HALVE ITSELF IN SUCCESION! 

CLT * 10# 

CL * CUT# 

T * X + H# X * FINAL* CONSTANTS! 

COMMENT MAIN KUTTA-MERSON STEP LOOP # 

LIFOR T*1 STEP H UNTIL FINAL 00 
BEGIN KM) F C T“H# Y# FO 3 # 

FORI Y 1 C I ] * F0tnxH3 + YU]# FTT-2xH3* Yl# Ft)# 

FORI Y 1 [ I ] <■ CFO[I]+Fim )xH 6+Y[I]# FCT-2*H3# Yl* FD! 

FORI Y1CI3 < (Fl[IIx3,0+FOEI])xH8+YCI3! FCT-H2# Yl# E23* 

FORI Y1EII * CF2inxA,0-FlC IIx3.0+F0[n )XH2 + Ytl3* FCT* Yl# FIT# 
FORI Y2CU * <F2[I3 xa,0+F1[ I3 + F0ni )XH6 + Y[I3! 

COMMENT DUES THE STEP SIZE H NEED TO BE CHANGED ! 

IF EPS^O THEN 
BEGIN OBL <• TRUE# 

FOHI BEGIN ERR*ABSCY1[I]-Y2[I3)xO, 2! TESUABSCYl 1 1 3 XEPS! 

IF LRR>TEST AND ERR>AB THEN COMMENT HALF H# 
bEGlN H «• H2# T4-T-H2! 

IF CCU*CU-1)<0 THEN ERROR * TRUE* 

IF STEPSIZE THEN WR I TE ( MSSG* H * T ) * 

IF T+H=T THEN BEGIN X*T* ERROR * TRUE! GO ID RETURN! 
END# 

CONSTANTS! GO TO KM! 

EnD# 

IF 64 * OxEHR>TEST THEN DBL * FALSE! 
tND* 

IF DBL AND H < HH THEN BEGIN H * 2.0xH! 


a Roll-Yaw 

Figure F 6. Solution of Optimal Linear Acquisition Equations m Backward 
Time 
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IF STEPSIZE THEN WRI TEC MSSG, H/ T ) ! 
t»U +■ CUT/ 

CONSTANTS END DOUBLE H! 

END J 

FORI YU] «• Y2CI3! 

END KUTTA MERSON LOOP! 

if eps=o then go to return; 

COMMENT NUW BE SURE TO HAVE T = FINAL ! 

HC <- H/ H «■ FIN AL"* ( T“H ) ! 

IF ABSlH)>AesCFINAL)xl,A55i9l5228^-U THEN 

BEGIN T f FINAL; EPS <• Q! CONSTANTS; GO TO L END! 

RETURNI END KUTTA flERSON! 

PROCEDURE DtlTAU/Y/DY]/ 

VALUE TAU/ REAL TAUj 
ARRAY Y/DYCO]/ 

BEGIN 

REAL A1,A2/A3/ A4/TAIJB; 

TAUB* TAU-TF-THETAO/ CT* COSETAUB)! ST «• SIN(TAUB); 

AH l + 2xE*CT/ A2<- 2x£xST/ A3<- 4+13x£xCT; A 4 <- l + 4xExcT/ 

COMMENT Yl = xi/ Y2=X?/ Y 3=X 3/ Y4=X4/ Y5 = P1/ Y6=P2» Y7=P3/ Y8=P4/ 

Yl«- YU]J Y2* Y C23 / Y3<- Y C 3 3 / Y4«- Y[43/ YS* Y [ 5 ] / Y6«- YE63! Y7<- YC73! 
Y8<- Y E 8 3 / 

Vl<- IF ABS(Y6)>1,0 THEN N1 xSIGN(Y6) ELSE 0; 

V2 *■ IF ABSIY83>1.0 THEN N2xSIGNCY8) ELSE 0; 

DYE 1 ]4-"Y2; 

D Y t 2 ]<■ KlxAAxYl-A2xY3“K3xAixY4»CxAtxEXP(KxEx(CT-i) )-Vi/ 

0 Y [ 3 3 «-"Y 4/ 

D Y [ 4 ]«• A2xY1+K4xA1xY2-K2xA3xy3“V2, 

DY[51<-“KlxA4XY6~A2xY8/ 

D Y [ 6 ] «- Y5-K4XA1XY8, 

D Y C 7 3 <- A2xY 6 +K2xA3xY8/ 

DYC83* K3XA1XY6 + Y7; 

QYE9 ]*• ABSC VD + ABSC V2)/ 

END DE> 

COMMENT INITIAL CFINAL) CONDITIONS j CLOCK! 

START. READCE/K1,K2/K3/K4/M/N2/X10/X20/X30/X40>P10/P20/P30/P4C/C/K, 
TF/THETAO)CFINISH] J 
T AU* TC03* Y E 9 ] <• 0/ 

Y E 1 3 <■ XI [U3<- Xio; 

Y E 2 3 *• X2C0J* X20/ 

YC33* X3 1 0 3 «■ X30! 

Y C 4 ]«- X4C0J«- X40/ 

y e 5 ]<■ pio; 

YE6K- P20! 

Y [ 7 3 «• P30/ 

Y C 8 3 <■ P40/ 

WRITE (LABL)/ 

WRITEERESL/TAU/X10/X20/X30/X40/P10/P20/P30/P40); 

COMMENT CALCUtATlNG WITH KUTTAMERSOn ANO LOADING ARRAYS! 

FOR L* 1 STEP 1 WHILE LS150 AND ABSC YE2 3 ) <1 , 850 AND ABSE YE 4 3 ><1 , 850 DO 
BEGIN 

KUTT AMERSUNC 9/ TAU/ 0 , 0 l / Y/ OEz £-4/ 8-5, NUTS/ FALSE)! 

IF NUTS THEN WR I TEC <"STEPS IZE WAS CUT AT LEAST TEN TIMES BUT KUTTAMERSO 
N C0NTlNUED/T="/F6,3/>/TAU); 

XI C L 3 «• Y E 1 3 / 
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X2ILI* YUJJ 
X3IL3* YC3J; 

X<l CU Y[4J, 

TU3 * L/ 

PRINT «■ If 

IF L MOD PRINT =o THEN WR ITE C RESL, T CL I > Y [ 1 1 » Y C 2 3 * Y l 3 1 > Y [ 4 I » Y [5 I * Y 1 6 ] > 

YC7], YC8] 

END CALCULATING AND LOADING LOOPJ 
J<- Y [ 9 ] j 

WR1TECCPAGE],PARM J .E,K1,K2,K3>K4,C, J)i 
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. . . . STANFORD B5500 ALGOL — 07/05/66 VERSION 201/67 

BEGIN COMMENT OPTIMAL PITCH ACOUISITION IN backwards time* 

SAVE 

FILE OUT PLOTTER 2U,765);REAL ARRAY A C 0 : 1022 3 , A2 C 0 t 2, 0 : 23 , C C 0 * 9 3 , S YMA C 0 
! 1123* SYmBI -15*63 ]>ABCD CO !2 3* REAL 0LTH, 

PROCEDURE PLOTCX*Y*IC);VALUE x,y,ic;real X»Y;iNTEGER ic;begin stream pro 
CEDURE WRlTCFlL,lFL»A);BEGIN SI*A;DI*IFL; 17(60(SI*SI+2,DS*6CHR) )J release 
(FIL 31END WRITJOWN INTEGER TT> I * NPX, NPY, B A* T, WlDX, PEN,* OWN BOOLEAN BQOL'I 
NTEGER J,K,DX,DY,IY,NR,NT,nC,II2,NA,LABEL FINSH, 0N'L1!F0RMAT TORMC 1022 A6 
);IF ABS ( IC )>3THEN BEGIN BA*0,WIDX*IF IC<-3THEN 2P2ELSE 1P2;PEN*2, A[03«-» 
444444";aC13*"444433",AC23*"333332",G0 TO FINSH,*END,IF A 8SCIO/PEN THEN 
Go 0NiPEN«-5-PENl Atn*-(IF BOOL THEN"66666l»ELSE TT+2A955)+PEN+PEN;B00L*TR 
UE, FOR K*9XPEN STEP-lUNTlL 1500 BEGIN A C 1*1 + 1 3 *"666666";iF I>1018THEN BE 
GIN AC l0i93+"340000",WRIT(PL0TTER>PL0TTERCO), AC 03 ), I*3lG0 ONI END, END, ON t 

NA*REAL<DX«—NPX + CNPX<-wiDXxX)>0)+REAL{DX>0 > * iy+Realcdy+-npy*cnpy*widxxy )* 

0) + REAL(DY>0 ), IF ABS( OX )>ABS C D Y )THEN BEGIN NR*ABS( DY > , NC+NT* ABSCDX ) M 12* 
A2[NA,1}END ELSE BEGIn NR* ABS CDX ) } NC*NT* ABS( DY ) i I I2+A2C 1, I Y 3ENDM Y* A2C N A 
,IY3,NA*nT DIV 2lNT*NT-NR,Ll*IF(NC<-NC-l3>0THEN BEGIN IF NA^NT THEN BEGIN 
T*IY,Na«-NA-NT;eND ELSE BEGIN T* 1 12 ,*N A*N A+NR ENOJIF BOOL*NOT BOOL THEN B 
EGlN ACI3*T+TT,IFCI*I+1)>1019THEN BEGIN AC 10193 *"340000"! WRIT C PLOTTER, PL 
OTTERCO), AC03 ), I*3,END>G0 TO L1,END,TT*"006006"&TC15'33;93;G0 to lijend; 
IF I C<OTHEN BEGIN IF BOOL THEN I*I-1ELSE AC 1 3 *TT + "000660"! A C I *1 + 1 3 *"3400 
00",WRITECPL0TTER,T0RM,F0R K*OSTEP 1UNTIL I DO AC K 3 ) , FI NSH * WR I TE C PLOTTER 
,TORM, "444444", "444433", "333331”, CCCBA*(BA+1)M00 100)M0D 103SCCBA DIV 10 
3 C24! 36s 123, "133333", "334444", FOR K+OSTEP 1UNTIL 70d0"444444") , BOOL*TRUE 

,npx*npy*o, i*3;endiend, 

PROCEDURE SYMBOLCXO, Y0,HGT, BCD, THETA, N 3, VALUE Xo, Y0,HGT, THETA, N, INTEGER 
NJREAL Xo,YO,HGT, theta; alpha ARRAY BCDCo3, BEGIN real BINX, AC,W,OSC, ainx, 
I,0STS,MnVl,HI,II,MoV2,WC,CC,Lp,XA,YA,XA6,YA6,0WN REAL CTH, STHIL ABEL LOA 
D8 , DEFINE A=SYMA«, B=SYMB#, I F T h ETA/OLT h THEN BEGIN OLTH*THETA; CTH*CQg( TH 
ETA*THETAx. 01 74532925 n;STH*5 INC THETA), END! H I *HGTx. 1428571 42857, XA*CTHxH 
1; YA*STHxHi;XA6*XAx6,0, YA6*YAX6.0; IF NcOTHEN’ BEGIN M0V1+IF Ns- 12THEN 3EL 
SE 2;xO*xO-MOVlx(xA+YA), Y0*Y0-M0V1 x(XA-YA3,BINX*N,G0 to loadriend; AC*wC* 
CC*0;W*BCDC03;WHILE AC*AC+1<N do BEGIN IF CC*CC+1 >7THEN BEGIN W*BCDCWC*W 
C+l3;CC*lEND;BINX*W,Cl2:63,W*0&WC12.18:3o3,LOADBtOSC*BCBINX3.C33i53,OSTS 
*aCAINX*bCBINX3,C39:933,LP*3)II*0;F0R I* 1STEP lUNTIL OSC DO BEGIN IF II* 
II + liSTHEN BEGIN II * 1 ; OSTS*A C A I NX* A I NX+ 1 3 END; MOV 1 *OSTS . C 6 t 3 3 > MO V2*GSTS . C 
9‘33;0 StS*0&0STSC6J12:36],IF M0V1=7tHEN LP*3ELSe BEGIN PLOTCXO+MOVlxXA-M 
0V2 xYA,Y0+M0V1xYA+MQV2xXA,lP),LP*2,END,END;X0*X0+XA6;Y0*Y0+YA6,END,EN0; 

procedure number cx,y,hgt,flt,theta,ni; value x, y,hgt,flt,thet a, n; integer 
n;real x,y,hgt,flt, theta, begin real frac, BOOLEAN b, label SWORDIREAL STRE 
am procedure cvcxo);value xo;begin si*loc xo;di*loc CV,DI*DI+2,DS*6DEC e 
ND, RE AL STREAM PROCEDURE LZROC VN, B ), VALUE VN,B, BEGIN LABEL M;SI*LOC VN,D 
I *LOC LZR0;DS*WDS,SI*SI-6;dI*DI-6;5(IF Sc=" 0"THEN BEGIN DS*L IT" ",Sl*SI + 
1 END ELSE JUMP OUT TO M),M.SI*LOC B,SKIP 47SB,IF SB THEN BEGIN DI*Ol-l,D 
S*LIT"""END;END L7R0 ;fRAC*IF NS0THEN.5ELSE IF N=1THEN,05ELSE IF N=2THEN, 
0Q5ELSE IF N=3THEN. 0005ELSE IF N=4THEN. 00005ELSE. 000005, B*FLT<-FRAC; IF F 
LT*ABSCFLT)+FRAC>100000THEN IF FLT210000000R B THEN BEGIN ABCDC 03*"*+*** 
*", GO TO SWORD END,ABCDC03*LZR0CCVCFRAC*ENTIER(rLT)),B);iF NSOTHEN SWoRD 
•SYMBOLCX,Y,HGT, ABCD, THETA, 6) ELSE BEGIN ABCDC 13*CVCENTIER( CFLT-FRAC )xiOO 
000))+", 00000", SYMBOLCX, Y,HGT,ABCD, THETA, IF N<5THEN N+7ELSE 12 1 .rwn itmd. 
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PROCEDURE AXISC X*Y* BCD* NC, SIZE* THETa* YM IN* DY), VALUE X, Y*NC* SIZE* THETA* YM 
IN* DY* REAL X*Y* SIZE*THETA*YMIN*DY* INTEGER NC * ALPHA ARRAY BCOCO)*BEGlN RE 
AL TH*CTH*STH* I*XB* YB, X A* YA* I CTH, ISTH* YM AX* ABSV.EXPPJ LABEL L50* TH«-THETAx 
0. 017455; ICTH<-CTH«-COS(TH)* ISTH«-STH*-SINCTH)*XB*X* YB*Y* IF THETA<OTHEN BEGI 
N ICTH < ‘"ICTH*ISTH«-~ISTH END* PLOTC XA*X- . lxl STH* Y A*Y + . lx ICTH , 3 3 i FOR I*1STE 
P 1UNTIL SIZE DO BEGIN PLOT C XB* YB* 2 ) } PLOT C XB*XB+CTH* YB*YB+STH* 2 ) * PLOT ( X A 
«-XA+CTH*YA«-YA+STH*2)*END*IF nc<othen go to l50*ymax*-ymin + sizexdy;if absv 
«-ABSCYMIn 5<ABS( YMAX )THEN ABSV«-ABSC YMAX ) * EXPP^O; I «-l* WH I LE ABSV>9999 . 999D0 
BEGIN ABSV«-ABSVx,i,l<.ix.i;EXPPfFXPP”lEND*WHILE ABSV<0,999D0 BEGIN ABSV«- 
ABSVxlO* 0*1^1 0*01 EXPP«-EXPP+1 END* OY«-DYxl I YMAX<-YHAXxI*XA<-XB-. 15xi STH-. 53 

xcTH*YA*YB+.15xICTH-.53xSTH*F0R I«-SIZE STEP-IUNTIL ODO BEGIN IF ABSCYMAX 
)>0,001THEN NUMBERCXa*YA*0. 1*YMaX*THETA*3)IYMAX«-YMAX“DY ;XA*Xa _ CTH*Ya*YA" 
STH END*I*IF EXPP=OTHEN NC ELSE NC + 7I Y A«-C ( S IZE«-S I ZEx , 5 ) - . 06 xl ) xSTH+ , 33X1 
CTH+Y;Xa«-(SIZE- .06xI)xCTH-. 33xISTH+X;SYMBOLCXA#YA* .1A*8CD*THETA*NC )* IF E 
XPP = OTHEN GO TO L 50 *I*-(I- 6 )x,i 2 *XA<-IxCTH + XA*YA«-IxSTH+YA;ABCD[o:i < -"x 10 tt 
iSYMBOLCXA* YA* .14*ABCD*THETA*6)*IF EXPP = lTHEN GO TO L50 ; X A*X A+ .25xCTH“. 0 
7xSTH*YA«-YA + .25xSTH + . 0 fxcTHJNUM 8 ERCXA*YA* ,07* EXPP* THET A* 0 > J L50 1 END* 

PROCEDURE LlNE(X*Y* N*K)* VALUE N*INTeG£R N*BOOLEAN K*ARRAY X*YCO];bEGIn I 
NTEGER I*I3*A*B*C*IF k THEM BEGIN A*-0* B*1 * C*N-1 * END ELSE BEGIN A«-N-lJB*“ 
1 JC«-0*END*I3«-3JF0R I*A STEP B UNTIL C DO BEGIN PLOT CX [ I ] * Y C I 3 , I 3 ) } 1 3«-2EN 

d*k*-not k*eno; 

PROCEDURE SCALE{A*N*K*L*YMIN*DY)I VALUE N,L*K* INTEGER N*K*REAL l*DY»YMIN* 
REAL ARRAY AtO)*BEGIN REAL YMAX*INTEGER I*NMX* YMIN«-YMAX«-ACK-n *NMK«-N-K*E 
OR I*K+K-1STEP K UNTIL NHK DO IP A C I ) > YM AX THEN YMAX*A 1 1 ]ELSE IF At I ]<YH 
IN THEN YMIN«-AtI]>DY«'CYMAX~YMIN)/L*EOR H-K-1STEP K UNTIL NMK DO Atl3<-(At 
I]-YMIN)/DY*END* 

PROCEDURE SETPLDTTER*BEGIN PLDT(0,0*0.0*-1019)* ABCDCO] «->*START "*ABCDC13<- 
"CALIBR"; A8CD[23*"ATION "* SYMBOL CO. 0* 5.0* 0.28* ABCD* -90,0* 18) I ABCDCO 3 «-”Y = 
0 "* SYMBOL (0.0* -0 . 14*0,28*ABCD* 0 .0* 3 ) * PLOTC 0 . 8*0 ,0*3 )* PLOTC 7, 0*0 ,0*23* 

ABCDC03*"SET X "* SYMBOL C6,7*0,5*0.28* ABCD* 90 , 0, 5 ) * PLOT C 7 ,0* 0 , 0* 3 ) * PLOT ( 7 
. 0*28, 0,2)1 PLOTC 1 .0,28.0,2)* ABCDCO ]<-"Y=28 ”* S YMBOLCO . 0, 27 . 84 * 0 ,28 * ARC D* 

0.0* 4)* PLOTC 10, 0*0.0* -3)* END* 


FILL C C * ) W I TH 0 CT040404 040404, 0CT04 0404040405 *0CT040404040 406, OC TO 404040 
40 407, 0 CTO 40 4040405 0 4, OcT040404040505*OcT040 40 40 40506, OcTO 40 404040507,0c 
TO 40 404040604* 0CT0404040 40605* FILL A2 CO** I WITH OCT50500*OCT506 00* DCTS070 
0,FILL A2C1,*]WITH OCT60500,OCT60600,OCT60700*FILL A2[2**]WITH 0CT705O0, 
0CT7O60O,0CT70700,FILL SYmaC*]WITh OCT 10 304 14637 1706 , OCT 1 1 0000 0000000 * OC 
T10302027160000*OCT40000144463717,OCT6050000000000,OCT01103041433414,OCT 
344546 37 170600 » 0 CT 4 30 3 37 30 20 4000 *OCTOOOOOOOOOOOOOO, OCTO 110 304 1433404* OCT 
7470000000000* OCTO 3 14 34 4 34 130 10, OCT 10 6 17 37460000, OCTO 60747 2 120 0000 *0CT34 
4 34 130 100 10 3* DC T 14 344546 37 1706* 0CT05 140000000000* OCTO 110304 146 37 17 *0CT60 
4 1333440000, OCT Hl5 1404 14 1202, 0CT42323 1 35344404 , QCT40 lOO 105 164655 , 0CT5 1 4 
2443525 14 12* OCT 213 14200000000, OCT iO H2120 107022* 0CT234546 37 170600, OCTlll 
22221 1170 14,0CT 15 25 24 1 40000 00, 0CT02 4 40600000000, OcTO 14 17006440 200, 0 CT2 12 
5230 34 30000, OCTOOO 3 4 346 37 1706, QCT34 34000000000, 0CTO407 3746 45 34 04, OCTO 304 
14 3340000* OCT 4 24 13° lOO l06i7*0CT 37464500000000* OCT0O07 37464 13000, DC T 470 70 
4 34040040 * OCT470704 3404 0000 *OCT 334 34 1 30 100 106, OCT 1737464500000 0* OCTOOO 70 
4it4474000*OCTl0302027173700*OCTooOlH10000000*OCT40202747000000*OCT362 71 
706O5403l*OcT4220i0Ol03 1400*0 cT 40 312329 364700, OcT 450 34 lOOOOOOOO* 0cT4 40 4 1 
5 1 3040000, DC TO 1 452 3054 1 0000, OCTO 1 10 304 1 470000, OCTOOO 70 34725400 0,OCT070 04 
000000000, OC T0007 2 3 4740 0000, OC T0007 4047000000, OCT 10 30 4 146 3770 36, OCT 477 0 3 
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7 1^0601 10 ,0CT0007 374645 3404 ,0CT 100 106 \7 37464 1, OCT 30 1070 22400000, OCTOOO 7 3 
7 4645 3404, OCT 344 34000000000, 0CT02 324 3 34 1405 16, 0CT462627 20000000, acTO 1452 
3054 1 2 30 3, OC T4 3 2325 2 10000 00, OCTO 34 300000 00000, OCTOO U23 24 160700, OCT 102 12 
212112170, 0CT2425 15142400 00, OCTO 14170 4 20 4 4600, 00700000000000000,00700470 

000000000, 007011030*14334 14, OCTOBO^ 17374 6 0000, OCT2027o747oOOOOO,OCTo7o 11 
030414700, OCT07204700000000,QCT0700 2440470000, OcTOO 4770 07400000, 0CT072 *4 
724200000,0CT0747O04070l434,0CT0OOC0000000000,0CTl02123l3l2220O,0CT00477 
0 160607 i7, OC T167041313040 4 l, 0CT0246 354505 35l 3, 0CT34 30000000000, OCTO 24270 
04*40000, QCT07272000000000, OCT 14167036340000, 0CT00404404000424 j 0CT220000 
00000000, OC T 004422 044 02200, OCT 00* 000440* 4422, 0CTO04 0044 4 002200* OCT 24 34 43 
41 30 100 1,0CT 31 42422000000, 0CT2**22002242 200, 0CT20224 42? 042200, 0CT22242 20 
0000000, 0CT240 14 124220000,OCT24024224202200,DCT4433 1304 131 100, 0CT1 131403 
1 3 32200, OCT 10 36500464 *070, OCT 22220000000 000, OCTO 24 22220 24 2200, OCTOO 440 44 
0220000, nCT 00 44 222*202204, oCT 40220242220 000, FILL SYmB[*)WITH 0CT30157, OC 
T 12156, OCT 14 155, 0CT22l53,0CT32l5l, OCT 14 150, OCT 12147, OCT 06146, OCT 141 45, DC 
T i4 144, 0cT26 1 42, OcT 14 1 41,0 C T 16140, 0cTl4l 37, 0CT2O1 35, 0cT22OO0,0CT 12002,0c 
T22003,0C7 32005,0CT 14007 , ocT 220 1 1, OCT300l3, OCT 120l5, 0CT400 16, OCT 30021, OC 
T34023,0cT42O25,0CT32030»0CT2o032,0cT06O34,0CTl4035,0CT12036,0CT24037,0C 
T 3004 1,0 cT 24O4 3, OcT 16045, OcT 16046, CcTi 4047, 0cT26050,0cT 14052, OcT 1405 3 , Oc 
T 1 2054, OCT 10055, OCT 32056, 0CT1 4060, CC TO 6061, OCT 12062, OCT 1206 3, OCT 1 2064, OC 
T 14065, OCT 06066, OCT 12067, OCT 10070, OCT 34071, OCT 1607 3, OCT 30074, OCT 24076, OC 
T26l0O,0CT26l02,0CTO4l04,0CTl4l05,0CT30106,0CTl*llO,0CTO0111,0CTO4112,0C 
T30113,OCT10115,OCT14H6,OCT06117*DCT12120»OCT12121,OCT12122,OCT16123,OC 
Tl4i25,0cT34l26,0cT22l3O,OcTl2l32,0CTlOl33,OcT12l34i0LTH«-^l0; 

PLOT C 0, 0, 1019); 

BEGIN 

REAL CT,ST,X,TAU,TF,THETAO,E,K3,N3,V3,X50,X60,P50,P60, J,C,K, 

Y1,Y2, Y3,Y4,X5MIN,X6MIN,0X5,DX6, 

INTEGER L, PRINT, 

ARRAY YCO»5],T,X5,X6CO:500), 

800LEAN NUTS) 

LABEL START, FINISH, 

ALPHA ARRAY HOZ, VERIO t 13 ; 

FORMAT LABLC"TAU",X3,"X5",X5,"X6",X5,’- , P5",X5,"P6"); 

FORMAT RESL(I3,4CX2,F5.2)), 

FORMAT PARMC"E=",F5.2,X2,"K3-",F6. 3, X2, "C = ", El 0 , 2, X2, " J=”, F5 . 2 ) ; 

procedure kuttamersoncn, x,hh, y, f, eps, ar , error, stepsize); 

VALUE N, HH, EPS,Ab , STEPSlZEJ INTEGER N; REAL X,HH, EPS,A8 ; 

R^AL ARRAY YCOl, PROCEDURE F* BOOLEAN ERROR, STePSiZe; 

COMMENT 1 VERSION OF 660518 660722 

EPS AND AB ARE THE RELATIVE AND ABSOLUTE ERROR BOUNDS RESP. 
STePSIZe TRUe TO write stepsize when changed 
STEPSIZE FALSE FOR NO OUTPUT 

ERROR Is SET TRUE IF STEPSIZE BECOMES TOO SMALL ELSE FALSE; 
begin comment kutta merson integrates a system of n first order 

ORDINARY DIFFERENTIAL EQUATIONS. SEE L. FOX, "NUMERICAL 
SOLUTION OF ORDINARY and PARTIAL differential 
EQUATIONS*', P. 24, PERMAGqN PRESS, 1962 / 

OWN REAL HC, FINAL, H2, H3, H6, H8, ERR, TEST, T, H; 

OWN INTEGER I,CU,CUT; OWN boolean DBLi label t, km, return; 

OWN real ARRAY Yl, Y2, FO, FI, F2t0j30]; 

comment except for hc, the own variables are for speed only; 

FORMAT MSSGC"THE step SIZE IS NOW", R12.5," AT Ts”,R12,5)J 
DEFINE FORI = FOR 1*1 STEP 1 UNTIL N DO #, 

CONSTANTS * H24-H/2.0, H3*H/3,0; H6*H76,0; h8*H/8.0 #; 
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comment check for initial entry and adjust h if necessary / 

ERROR «• FALSE/ 

H *• HH / 

IF N=0 THEN BEGIN HC,«- H/ GO TO RETURN END! 

IF H=0 THEN GO TO RETURN/ FINAL «- X+H) 

if hc=o Then hc *■ h, 

IF EPSXO AND ABSCH)>ABS(HC) THEN 

IF SIGNCH)XSIGN(HC) THEN H <• HC «• -HC ELSE H «• HC/ 

COMMENT: CUT IS THE NUMBER OF TIMES THAT THE STEP SIZE IS ALLOWED TO 

halve itself in succesion; 

CUT 10 ; 
cu «• cut; 

T «• X+Hi X <- FINAL/ constants; 

COMMENT MAIN KUTTA-MeRSON step loop ; 

L» FOR T«-T STEP H UNTIL FINAL DO 
BEGIN KM 5 FCT-H/Y/FO); 

FORI Yim * FOmxH3+YCI3; FCT-2XH3/ Yl/ FI)/ 

FORI Y 1 C 1 3 «■ (FOC I]+Fim )xH6+Yt I]/ F(T-2 xH 3/ Yl, F1)J 

FORI YlII) *• (Fin)x3.0 + F0Cn)xHB+YCn/ FCT-H2/ Yi/ F2); 

FORI YlII] «- (F2Cl]x4.0-FlII]x3.0+F0Cn)xH2+YCn/ FIT, Yl/ Fl)/ 
FORI Y2CI] *■ (F2Cnx«.0 + Fltn+F0[I])xH6+YCn/ 

COMMENT DDES THE STEP SIZE H NEED TO BE CHANGED i 
if eps/o then 
BEGIN DBL «- TRUE/ 

FOrI BEGIN ERR«-ABSCY1CI3-Y2EI3)xO,2/ TeST*A 8S( Yl [ 13 ) xEPS; 

IF ERR>TEST AND ERR>AB THEN COMMENT HALF H; 
begin h «■ H2; t«-t-h2/ 

if ccu*-cu-n<o then error «■ true; 

IF StEPSlZE THEN WRITECMSSG/H/T)/ 

IF T+H=T then BEGIN X«-T/ ERROR * TRUE/ GO TO RETURN/ 
END/ 

CONSTANTS/ GO TO KMI 
END/ 

IF 64.0xeRR>TEST THEN DBL «• EALSE/ 

end; 

IF DBL AND H < HH THEN BEGIN H *■ 2.0xh, 

IF STEPSIZE THEN WR I TE( MSSG, H/ T 3 / 

cu «■ cut; 

CONSTANTS END DOU8LE H; 

end; 

FORI YCI3 «- Y2II3/ 

END KUTTA MERSON LOOP; 

IF EPS=0 THEN GO TO RETURN/ 

COMMENT NOW BE SURE TO HAVE T = FINAL ) 

HC «• H/ H «■ FINAL-CT-H)/ 

IF ABSCH)>ABSCFINAL)xl.45Sl9l5228P-li THEN 

begin t «■ final; eps <■ o, constants; go to l end; 

RETURNS end KUTTA MERSON/ 

PROCEDURE CLOCK;BEGIN OWN INTEGER TEMPUS, TEMPUS 1 / 

FORMAT FMTKX97/ "DATE: ", A2,A2,A2), FMT2<6 C 

“ - - ")/ "ELAPSED TIME WAS"/ F7,3# " SECONDS - - - - " 

" total time was"/ f8.3, » seconds,")/ 

IF TEMPUS X o THEN WRITECFMT2/ “(TEMPUS -(TEMPUS *■ TIMEC2) ) ) / 60.0 

/-(TEMPUSl - TEMPUS) / 60,0) 

EL s E BEGIN WRITECFMTl,CTEMPUS«-TlMEC5)>.I36»123, TEMPOS. 124:123, 
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TEMPUS »Cl2tl23!)* TEMPuS«-TEMPUS 1 <-T IME! 2 3 3 ENDiEND CLOCK/ 

PROCEDURE 0E!TAU*Y,DY)/ 

value tauj real tau; 
array Y,DY!03* 
begin 

REAL A3* T AUB/ 

TAUB* TAU-TE-THETAOI CT«- COS(TAUB); ST* SlNCTAUB)* A3*1+3 xExCT, Y1*YU3* 
Y2*Y[23;Y3*YI33*Y4*Y[43* V3*IF ABS(Y4)>1,0 THEN N3xSlGN<Y4) ELSE 0, 

DYC 1 ]*— Y2* 

DY C 2] *• 3xK3xA3xY1+2xexST+CxC 1+2xExcT )XEXP!KxEx!CT-1 ) )-V3; 

DYI 33*-3xK3xA3xY4; 

DY C 4 3 «• Y3J 
DYC5]* ABS(V3)*' 

END DEi 

COMMENT INITIAL (FINAL) CONDITIONS; clock; 

START* READCE*K3*N3*X50*X60*P50*P60*C*K*TF*THETAO)tFINISH3; 
TAU*T!0]«-YC53*0* 

Yin* xbcoi * xso; 

YC23* XAC03* X60; 

Y E 3 3 4- P50* 

YC43* P60; 

WRITE CLABL3* 

WRITE !RESL/TAU*X50*X60*P50*P60)* 

COMMENT CALCULATING WITH KUTTAMERSON AND LOADING ARRAYS! 

FOR L* 1 STEP 1 WHILE L<150 AND ABS C Y C 2 3 )<1 . 850 DO 

begin 

KUTT AMERSGN { 5* TAU* 0,0 i*Y*DE*P"4* P-5* NUTS* FALSE); 

IF nUTS THEN WRITECchSTEPSIZE WAS CUT AT LEAST TEN TIMES BUT KUTTAMERSO 
N C0NTlNUED*T= n /F6,3/>/TAU)/ 

XS CL3 * YC13; 

X6CL3<- YC23! 

TCL3 ^ l) 

PRINT <*1* 

IF L MOD PRINT =0 THEN WRITE! RESL* TIL) * Y I 1 3 * Y E 2 ) * Y I 3] * Y C 4 3 3 i 
END CALCULATING AND LOADING LOOp; 

J «• Y C 5] * 

WRITE! IPAGE3*PARM*E*K3*C* J); 

COMMENT PLOTTING X6 VS X5 * CLOCK* 

XSILI^XAILI^-R* X5!L+l3«-X6CL+i3«-2; 

SCALECX5/L+2/ I* 4* X5MIN* DX5 ) * 

SCALE!X6>L+2/l/4*X6MlN*DX6); 

HQZI03* " X5 "* 

VERC03* " X6 «; 

PL0T(0*l,-3)j 

AXIS!0*0*H0Z*4*4*0»X5MIN*DX5); 

AXIS!0*0*VER*4*4,90*x6MIN*Dx6)I 

LINE!X5*X6*L/TRUE)* 

PLOT CIO* -l* -3 3* CLOCK! 

COMMENT IF MORE DATA CARDS ARE TO BE READ THEN; 

GO TO start; 
finish; 
end; 

PLOTcO/0/1019); 

END, 
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